format: text

text:
# STJ Observations of the Eclipsing Polar HU Aqr

## 1 Introduction

HU Aqr is a member of the AM Her sub-class of cataclysmic variables, also called polars because of their highly polarised optical emission. Polars are characterised by a white dwarf primary with a strong magnetic field ($`10240`$ MG) and no accretion disc (see Cropper 1990 and Warner 1995 for reviews). The secondary is a main sequence star which fills its Roche lobe, resulting in an overflow of material from the inner Lagrangian point to form an accretion stream. The stream is initially on a ballistic freefall trajectory, until at some radius the magnetic pressure of the white dwarf magnetic field is approximately equal to the ram pressure of the stream. At this point the stream is threaded onto the field lines of the white dwarf and accretes onto a small area near one or both of the magnetic poles.

HU Aqr is an eclipsing polar with a period of $``$ $`125`$ mins, lying just below the period gap (see e.g. Howell, Nelson & Rappaport 2001). Harrop-Allin, Hakala & Cropper (1999b) used the eclipse of the white dwarf primary by the secondary to establish the distribution of brightness along the stream given a number of parameters defining both the geometry of the system and the physical parameters of the model. The stream emission was calculated using model light curves which are optimized by a genetic algorithm, based on a method first employed on photometric eclipse profiles of HU Aqr by Hakala (1995). The technique was developed and improved by Harrop-Allin et al. (1999a) and applied to observations of HU Aqr in both a high and low accretion state (high- and low-mass transfer; Harrop-Allin et al. 1999b and Harrop-Allin, Potter & Cropper 2001).

This method of eclipse mapping has also been applied to emission-lines (Sohl & Wynn 1999; Vrielmann & Schwope 2001; Kube, Gänsicke and Beurmann 2000). Kube et al. (2000) used HST FOS spectra of (another polar) UZ For to map the accretion stream line emission at Civ $`\lambda `$ 1550Å, and Vrielmann & Schwope (2001) used the emission line light curves of H$`\beta `$, H$`\gamma `$ and He ii $`\lambda `$4686Å of HU Aqr. Both these line-emission methods use a three-dimensional stream so data from the entire orbit of the system is fitted and the model can reproduce the out-of-eclipse features, such as the pre-eclipse dips (Watson 1995). They can also reproduce the effect of differing brightness distributions between the irradiated and un-irradiated faces of the stream, and this effect is not reproduced in the model of Harrop-Allin et al. (1999a and 1999b). However, the Harrop-Allin et al. technique has the benefit of being sensitive to the total emission from the stream (not just the line emission). The technique also has the advantage of access to high signal-to-noise ratio data, which is crucial in producing robust fits.

The light curves presented here were obtained using the S-Cam2 Superconducting Tunnel Junction (STJ) camera, developed by the ESA Astrophysics Division at ESTEC. The camera is the second prototype of a new generation of detectors that record the energy as well as the position and time of arrival (to within $`5\mu `$s for this particular detector) of the incident photons (Rando et al. 2000). The application of STJs to optical photon counting was first proposed by Perryman, Foden & Peacock (1993), and has since been applied to observations of the Crab Pulsar (Perryman et al. 1999), the magnetic cataclysmic variable UZ For (Perryman et al. 2001) and to quasar spectroscopy (de Bruijne et al. 2002).

The detector itself consists of a liquid helium cooled array of $`6\times 6`$ pixels, each being $`25\times 25`$ $`\mu `$m2, corresponding to $`0.6\times 0.6`$ arcsec2 per pixel and a field-of-view of about $`4\times 4`$ arcsec2. The pixels are sandwiches of superconducting tantalum with a thin insulating layer in the middle and the whole device is cooled well below the superconductor’s critical temperature (about 0.1Tc). An incident photon then perturbs the device equilibrium and as the energy gap between the ground state and the excited state is only a few meV, a large number of free electrons is created by each photon, this number being proportional to the photon energy. This is in contrast to normal optical CCD semiconductors where the band gap is $`1`$ eV, and photon absorption results in typically only one free electron being created.

The S-Cam2 instrument is particularly suited to observations of cataclysmic variables. The high-time resolution and simultaneous observations of spectral and intensity variations is ideal for eclipsing systems with orbital periods of the order of those in cataclysmic variables. The intensity variations over the rapid ingress and egress can be probed directly, and the rapid variations of the intensity of the accretion stream can be used to provide information on the possible mechanisms of emission along the stream.

## 2 Observations

S-Cam2 was mounted at the Nasmyth focus of the William Herschel Telescope during 2000 October 2/3 and 3/4. Table 1 gives the cycle numbers (relative to the ephemeris of Schwope et al. 2001), start time in UTC, time of mid egress and the observation length in seconds.

A total of five eclipses of HU Aqr were recorded. The ingress was missed for eclipse 29982 (ephemeris of Schwope et al. 2001), and only just recorded for 29993. The relatively high count rate of HU Aqr caused some difficulties for the data acquisition system, limiting the duration of the runs. The first three eclipses exceeded the data acquisition system limits, resulting in loss of absolute time reference. Eclipse 29994 was therefore taken with a neutral density filter with a throughput reduction factor of 10. This is the most complete of the eclipses, but suffers from reduced signal-to-noise ratio. The data gap in cycle 29982 was caused by the problems with the data acquisition system noted above. The seeing was in the range 1 to 1.5 arcsec, except at the beginning of the run for cycle 29982 and throughout cycle 29983, when it was 2 arcsec, which led to some spillage of light from the $`6\times 6`$ array. We have therefore not used those data.

Observations of the spectrophotometric flux standard BD+28 4211 were taken to calibrate the data. They were made through a neutral density filter with attenuation factor 100 and also suffered from poor seeing.

## 3 Data reduction

The data reduction process employs specific pipeline processing of the S-Cam2 data (de Bruijne et al. 2001), based on the FTOOLS suite of software (Blackburn 1995) and a full description is given by Perryman et al. (2001) for their observations of UZ For.

Once reduced, the data can be split a posteriori into different energy bands (wavelength ranges). The intrinsic STJ resolution is such that $`\lambda `$/$`\mathrm{\Delta }\lambda 9`$, but here we use only three bands to represent ‘red’, ‘yellow’ and ‘blue’ light. The energy bands are selected so that roughly equal numbers of events occur in each. The wavelength range is limited by the atmosphere at short wavelengths and the use of optics to suppress infrared photons at long wavelengths, and corresponds to approximately 340$``$680 nm. The wavelength ranges for the blue, yellow and red bands here are 340$``$470 nm, 470$``$550 nm and 550$``$680 nm.

As part of the pipeline processing, the data were flat-field corrected using a single map derived from sky observations. The data were then corrected for atmospheric extinction using the standard La Palma tables of extinction values as a function of wavelength and air mass, and finally sky-background subtracted. The pipeline processing is applied to a light curve in arbitrary time bins of 1 s. The extraction of the background-subtracted light curves used the entire $`6\times 6`$ array for the object and the two corner pixels (1,1) and (6,6) for the background. The background light level is then taken as the mean level during eclipse and this is subtracted from the source light curve.

Figure 1 shows the white light curves for eclipses 29982, 29994 and 29995 before subtraction of the background, together with the backgrounds taken from the two corner pixels (1,1) and (6,6).

Once the data had been calibrated and reduced, we folded the data on the orbital period. The linear ephemeris of Schwope et al. (2001) was used, which defines inferior conjunction of the secondary as $`\varphi =1.0`$. The UTC times are transformed to TDB (at the solar system barycentre, i.e. including light travel times). For cycles 29982 and 29993, where there is no absolute time reference, the egress is aligned to be at the same phase as in cycles 29994 and 29995.

Figures 2 and 3 show the white, red, yellow and blue light curves, with the colour ratios yellow/red and blue/yellow, for cycles 29993, 29994 and 29995. The secondary star has been subtracted from the colour ratios (but not the light curves) so that the change in the ratios is due wholly to the stream, after accretion region ingress.

We extracted intervals of good seeing from the standard star observation to determine the countrate in the yellow band, corresponding most closely to the $`V`$ band. This gave a zero point magnitude for 1 count/s in this band of 24.0, taking into account the (assumed) factor of 100 from the neutral density filter. The corresponding maximum brightness at $`\varphi =1.1`$ in Figure 3 is $`V=14.7`$. This indicates that HU Aqr was in a high accretion state at the time of these observations (see Schwope et al. 2001).

## 4 The Light Curves

The light curves in Figures 2 and 3 show a number of features that are characteristic of an eclipsing polar system. The most prominent feature is the eclipse itself, which starts with the limb of the secondary star eclipsing the accretion stream. At about $`\varphi `$=0.964 the bright accretion region on the white dwarf is eclipsed in a few seconds (the white dwarf is also eclipsed around these phases, but is much fainter). After this only the secondary star is visible, and the sequence is then approximately reversed on egress.

### 4.1 Pre-eclipse dip

Prior to the eclipse of the white dwarf, we observe a pre-eclipse dip. Figure 3 shows this dip centred at $`\varphi 0.872`$. The dip is caused by the eclipse of the accretion region on the white dwarf by the magnetically confined section of the accretion stream. Consequently the phase of the centre of the dip is directly related to the azimuth of the coupling region, where the stream is threaded onto the field lines of the white dwarf. The centre of the dip in eclipse 29994 corresponds to an azimuthal angle of 46 which is within the range of values found by Harrop-Allin et al. (1999b) of 42 to 48 for their light curves.

Schwope et al. (2001) investigated the correlation between the phase of the dip and the soft X-ray luminosity, and concluded that the dip ingress occurs earlier in phase when the system is brighter. The ingress of the dip at $`\varphi 0.856`$ indicates that the soft X-ray luminosity is at the higher end of the range, as expected for the bright state of the system.

The dip is visible only in cycle 29994, the other observations being too short to cover the relevant phases. We cannot therefore relate its movement between successive eclipses to the change in stream eclipse profile as discussed in Section 4.3.

### 4.2 Accretion region eclipse

The soft X-ray data in Schwope et al. (2001) suggest that there is only one accreting pole. It is possible that there is a second accreting pole, as discussed in Harrop-Allin et al. (1999b). However, there is no evidence for such a pole in our eclipse profiles, so for the purpose of our modelling we assume accretion onto one pole.

Around $`\varphi 0.964`$ the accretion region has been completely eclipsed and the accretion stream is thereafter the dominant source of the observed brightness with a small contribution in the red band from the secondary. We show in Figure 4 the accretion spot profile on ingress and egress. The profile is constructed by subtracting successive intensities over 1 s time intervals (0.00013 phase) and so shows the rate of change of the light curves intensity. The profile on egress is the mean of cycles 29982, 29993 and 29995 and is therefore more reliable than the ingress profile, which uses only cycle 29995. The profile may be asymmetrical, with the leading edge (corresponding to later phases) of the cyclotron emitting region being brighter than the trailing edge. The accretion spot ingresses last for $`7`$ s and the egresses are $`8`$ s. This is unchanged from the durations measured by Case (1996) in the data used by Harrop-Allin et al. (1999b) and should be compared to the duration of $`1.3`$ s for the soft X-ray eclipses (Schwope et al. 2001), so that the extent of the cyclotron emitting region is $`15^{}`$ in the optical compared to $`3^{}`$ in soft X-rays. From soft X-ray absorption modelling, Schwope et al. (2001) suggest that the threading region is extended in azimuth, but the soft X-ray emission originates over only a small range of azimuths: the larger region seen in the optical is therefore likely to be the result of this extension in azimuthal threading. The cyclotron emitting region is significantly more extended than that measured for UZ For, where Perryman et al. (2001) found both regions extended over $`3^{}\pm 1^{}`$.

In those cycles where we have absolute time information (29994 and 29995) the egress takes place at $`\varphi =1.0413`$ in both cycles, 10 s later than predicted by the linear ephemeris. This indicates that the observed-minus-calculated (O-C) residuals in Schwope et al. (2001) Figure 4 are continuing to increase. In case there are systematic differences between the optical and soft X-ray egress times, we rephased the Harrop-Allin (1999) data in their Figure 1 on the new ephemeris, including the UTC-TDB time correction of 59 s appropriate to the Schwope et al. (1997) ephemeris epoch. The data were taken on cycles $`1322`$ in Figure 4 of Schwope et al. (2001), and the O-C of $`1.6`$ s is in agreement to within $`1`$ s with those determined from the soft X-ray data. The continuing increase in O-C residuals is an indication that the period in the linear ephemeris of Schwope et al. (2001) is slightly too short, and is in the opposite sense to that expected from their quadratic ephemeris which predicts a decrease in the O-C residuals. A period of 0.0868204212 days would eliminate the 10 s residuals we observe.

The duration of the eclipse from the centre of the accretion region ingress to the centre of its egress is $`\mathrm{\Delta }\varphi =0.0778\pm 0.0002`$ in cycles 29993, 29994 and 29995. This is the same as that found by Harrop-Allin et al. (1999b) of $`\mathrm{\Delta }\varphi =0.0779\pm 0.0002`$, and Schwope et al. (1997) of $`\mathrm{\Delta }\varphi =0.0779`$. The lack of variation in the duration of the eclipse is remarkable and implies that the accretion spot is at the same latitude for these observations. On the other hand, this is shorter by 6 s than the $`0.07820.0792`$ modelled by Schwope et al. (2001) from the soft X-ray eclipses, where the ingress was uncertain due to the low countrate and absorption.

For a particular mass ratio, the width of the eclipse determines the inclination of the system. Using $`q=M_2/M_1=0.25`$ (Schwope et al. 2001), the inclination is $`85.0^{}`$ (assuming $`M_1=0.9M_{}`$ and $`\beta =40^{}`$, $`\zeta =50^{}`$, where $`M_1`$ is the mass of the white dwarf and $`\beta ,\zeta `$ are the magnetic longitude and colatitude of the accretion region. These correct for the location of the accretion spot on the surface of the primary; see below).

### 4.3 Stream eclipse variations

Complete ingress of the stream occurs at $`\varphi =0.982,0.979,0.977`$ and $`0.982`$ for eclipses 29982, 29993, 29994 and 29995 respectively. During the period of complete eclipse the secondary is the only contributor to the light curve and provides a constant contribution which is greatest in the red band, as expected, with no contribution in the blue band.

Figure 5 shows more clearly that the observations reach total stream eclipse at different phases even though the spot ingresses/egresses both occur at the same time. This rapid change in the accretion flow$``$magnetic field interaction between successive eclipses was first observed by Glenn et al. (1994) who found a difference in the time it took for the stream to be completely eclipsed of more than a minute between their two successive eclipses.

The shape of the accretion stream eclipse profile in the observed light curves gives a qualitative prediction of where brightness enhancements can be expected in the accretion stream. When the stream is brighter at later phases in the stream eclipse, as in cycles 29982 and 29995 (Figure 5), there must be more emission in the threading region, as this is the only part of the stream still visible. For the same reason, if the final stage of the stream eclipse takes place at later phases, the threading region must have moved: thus in cycles 29982 and 29995, either the accretion stream is brought further out from the line of centres during the ballistic part of the trajectory, or the magnetically channeled stream rises further out of the orbital plane (we address this in Section 5 below). This could be caused either by a small change in the location of the accretion spot on the primary surface, moving the accretion spot further from the line of centres, or perhaps a change in the velocity of the stream as it leaves the L1 point. Likewise, if the accretion stream is fainter near the white dwarf, the system brightness will be fainter immediately after accretion region ingress. This may be the case for cycle 29995 (Figure 5).

### 4.4 Colour variations

For cycle 29995 there is a rise in the yellow/red ratio during the stream ingress (Figure 2), so that successively more red light is blocked from the stream as the eclipse progresses. This implies that the last part of the stream to be eclipsed, near the threading region, emits relatively less at longer wavelengths than the rest of the stream. This rise is not evident in cycle 29993, but maybe present at a reduced level in cycle 29994 (Figure 3). This trend is consistent with a threading region which is increasingly hotter in the three successive eclipses. At the same time as the large yellow/red ratio change in cycle 29995 there is no change in the blue/yellow ratio.

After eclipse when the stream close to the primary is uncovered first, the different yellow/red ratios in cycles 29993 and 29995 are consistent with a hotter stream near the white dwarf in eclipse 29993 than in 29995. An increase in the blue/yellow ratio after eclipse is seen in all three cycles (29993, 29994 and 29995) and also implies that the magnetically confined stream material is becoming hotter towards the threading region.

The pre-eclipse dip around $`\varphi 0.87`$ in cycle 29994 (Figure 3) is also slightly bluer than at earlier or later phases. This was also the case in the light curves of Harrop-Allin et al. (1999b).

## 5 Indirect imaging

### 5.1 Methodology

We have inferred from a qualitative analysis of the light curves that the trajectory and/or intensity of the stream varies over successive orbital periods, and we can identify brightness enhancements along the stream length. However, it is more difficult to infer the location of the brightness enhancements without a model for the stream trajectory. In particular, at stream ingress there is no information on which parts of the stream are eclipsed at which phases. We can use an eclipse mapping method to create a model accretion stream by placing model stream points along a predetermined trajectory, consisting of a ballistic part from the L1 point, coupling to a magnetically confined section where the material is threaded by the field lines of the white dwarf. The distance from the white dwarf where this transition occurs is the threading radius $`R_\mu `$, and is a user input to the model. The model system is rotated while the Roche lobe filling secondary eclipses the stream and white dwarf, with the brightness of the visible model stream points at each phase summed to form a model light curve. The resultant light curve is optimized with the genetic algorithm (GA) evolving the best fit light curve.

The ‘goodness of fit’ of a model light curve is measured with a ‘fitness function’ (see Harrop-Allin et al. 1999a for a full description). The GA adjusts the brightness points along the stream in an attempt to minimise this function, which consists of a $`\chi ^2`$ term and a maximum entropy term, which ensures the problem is not under-constrained. After the final stages of the GA, a more conventional line-minimisation routine (Powell’s method) is used to reach the final minimum of the solution.

The accretion stream itself is taken to be physically thin in that it does not eclipse the primary, and so features such as the pre-eclipse dip (e.g. Watson 1995) are not reproduced by the model. However, the brightness contribution of each stream point is taken as the sine of the angle between the line of sight and the tangent to the stream at that point, essentially a projection effect, thus optically thick. We also exclude the white dwarf from the model, although the egress of the white dwarf can be seen in the light curves immediately before the egress of the accretion region.

The exact method we apply for the modelling is different to that of Harrop-Allin et al. (1999b, 2001) in that we are not applying the technique to the complete eclipse light curves, nor can we interpret the results in the same manner. There are a number of reasons for this. Firstly, and importantly, our observations are all truncated, either in the ingress or the egress of the accretion stream. For observations where the egress is truncated (cycle 29995), the model has difficulty breaking the ambiguity which exists when assigning brightness to points eclipsed in the same phase interval, i.e. the brightness can be either in the magnetically confined region or the ballistic. For observations where the pre-eclipse light curve is missing, we lack information on the accretion region, and parts of the stream towards the secondary (cycle 29993). A further difficulty lies in determining the nature of the variable emission from the accretion region. This is particularly important at egress, where the brightness of the spot can greatly affect the brightness distribution towards the white dwarf. After the egress of the accretion region the light curve consists of a variable component, plus the stream, which changes rapidly over the course of the observation. We have therefore chosen to use a restricted modelling technique, which relies on the ingress of the accretion stream alone. This means that we have truncated the light curves, removing all phases up to and including the ingress of the accretion region, and those after $`\varphi =1.0`$. This is the only part of the light curve unaffected by the variable emission from the accretion region. We still have the problem of the ambiguity of the stream points, but we can make some progress in interpreting our results by referring back to the light curves and the colour ratios.

### 5.2 Fixed parameters and $`R_\mu `$

A number of parameters is required to produce a model light curve and reproduce the stream brightness distribution. These fall into two types: physical parameters such as the masses of the two component stars, and geometric parameters such as the location of the accretion region on the white dwarf (see Harrop-Allin et al. 1999a for details). The model is particularly sensitive to the exact value chosen for the parameter $`R_\mu `$, but we can constrain the value using the model stream geometry. If the value is too large then the model needs to assign a large amount of brightness to a few points. Conversely if $`R_\mu `$ is too small then no emission is assigned by the model to the points in the threading region. From brightness maps for a range of values of $`R_\mu `$ we can determine the end point of the ballistic trajectory and so provide a constraint on the value of $`R_\mu `$. We have found from our model fits of the light curves used here that the technique depends heavily on the data being of a sufficiently high signal-to-noise ratio, which is important because of the sensitivity to the value of $`R_\mu `$. Therefore, in order to determine the best fit value to use we further restrict the application of our model technique to cycles 29993 and 29995 (Figure 2).

Figure 6 shows model fits to the ‘blue’ light curve of cycle 29995 for different values of $`R_\mu `$. The stream maps highlight the dependence of the fits on a correct value for $`R_\mu `$. For a value of $`R_\mu =0.18a`$ the brightness of the threading region is low. This lack of emission, or ‘hole’, is caused by $`R_\mu `$ being too small - the data are incompatible with emission at the end part of the resulting long ballistic stream. On the other hand, if $`R_\mu `$ is too large, so that the model ballistic stream is shorter than it is in reality, a pile up of excess brightness at the threading region is seen, as in the stream map for $`R_\mu =0.26`$. The inappropriateness of the latter can be deduced from the poor fit to the region at $`\varphi =0.98`$ (corresponding to the threading region) for values of $`R_\mu `$ which are too large (insets to Figure 6).

For all our modelling we used fixed values of $`M_1=0.9M_{}`$, $`q=0.25`$ and $`i=85.0^{}`$ (Harrop-Allin et al. 1999b, Schwope et al. 2001). Table 2 gives the values of $`R_\mu `$ and the field orientation parameters for each cycle. We estimate the range on the values of $`\beta `$ and $`\zeta `$ as $`\pm 10^{}(1\sigma )`$.

The values of the three parameters in Table 2 are different from those found by Harrop-Allin et al. (1999b). We find $`R_\mu =0.220.26a`$ ($`1.3\times 10^{10}`$ cm to $`1.6\times 10^{10}`$ cm), whereas Harrop-Allin et al. (1999b) found $`R_\mu =0.170.20a`$. We can attribute this difference to an improved stream trajectory compared to that of Harrop-Allin et al. (1999a). This arises from an incomplete treatment of the appropriate forces in the frame of reference chosen by Harrop-Allin et al. (1999a). The difference is small, but for the angles used here the resulting $`R_\mu `$ is significantly different, with the Harrop-Allin et al. trajectory underestimating the values for $`R_\mu `$.

The decrease in the value of $`R_\mu `$ from cycle 29993 to 29995 implies that the stream penetrates further into the magnetosphere before threading onto the magnetic field lines (Figure 8). The effect of the change in the threading radius can be seen in Figure 5. The values of $`\beta 60^{}`$, $`\zeta 40^{}`$ found here are also larger than those used by Harrop-Allin et al. (1999b), but consistent within the uncertainties. Any change in these parameters would imply that the location of the accretion spot has moved further from the line of centres and the altered geometry of the magnetic field carries the accretion stream further round the white dwarf from the line of centres (see Section 5). They are in agreement with those found by Schwope et al. (2001) (their figure 6).

### 5.3 Model results

With the values of $`R_\mu `$ determined above, we can use the model to find the best fit brightness map, from which we can then determine the brightness per unit stream length along the stream trajectory. This will show the brightness independent of the actual length of the stream eclipsed in any particular phase interval. This is important as there are different lengths of ballistic and magnetically confined stream eclipsed in a given phase interval.

Figure 7 presents the results of the modelling in two columns representing the two cycles presented here and the three rows the three energy bands: red, yellow and blue. We have divided the stream into three sections, as illustrated in Figure 8, which enclose the section nearest the white dwarf, the threading region and the stream in between the two. The ordinate then represents the total brightness of points eclipsed in these intervals per unit stream length, that is points in the ballistic and magnetically confined sections of the stream eclipsed in that phase range. The results have been calibrated into energy units using the standard star observation in order to facilitate comparison between the different bands. However the calibration is only approximate due to a lack of a standard star on the second night of observations.

We can see from Figure 7 that there is enhanced brightness in the section containing the threading region in all bands in cycle 29993, but only in the yellow and blue in cycle 29995. The section containing those parts of the stream nearest the white dwarf are also bright in cycle 29993. However, in cycle 29995 this is only seen in the blue. The brightness of the different sections is consistent with the colour ratios from Section 4.4. In Figure 7 the threading region in cycle 29995 is hotter than that in cycle 29993 as the ratio of the yellow to red is greater for cycle 29995. This is again consistent with the colour ratios (Section 4.4). However, the blue to yellow ratio is similar in the two eclipse ingresses suggesting that the increase has occurred in both blue and yellow bands.

## 6 Discussion

### 6.1 Comparison with earlier results

The model fits of Harrop-Allin et al. (1999b) and Harrop-Allin (1999) for one pole accretion, show a general enhancement in the threading region and towards the white dwarf. Some eclipses (for example cycle 3723 and 3724 in Harrop-Allin 1999) show significantly greater brightness in the threading region of the white dwarf than others (such as the immediately preceding cycle 3722).

Enhanced brightness regions are also found by Vrielmann & Schwope (2001) and Kube et al. (2000) who apply similar modelling techniques. Vrielmann & Schwope (2001) apply their Accretion Stream Mapping technique to H$`\beta `$, H$`\gamma `$ and He ii $`\lambda `$ 4686Å observations of HU Aqr. They find a brightness enhancement at the threading region, but not towards the white dwarf, similar to cycle 29995, so implying that this may be absent in their emission lines. Kube et al. (2000) used the Civ $`\lambda `$ 1550Å line emission from UZ For, with a 3-dimensional stream model. They found three regions of enhanced brightness: one on the ballistic accretion flow, and two on the magnetically confined section. They suggest that the enhancements on the magnetically confined section are caused by irradiation of denser sections of the stream with a large area near the accretion region, and a smaller region near to the threading region. This may be the case in cycle 29993, where we find a bright stream near the white dwarf, and again near the threading region. Although Kube et al. find no enhancement actually at the threading region, this may be a result of an increase in the density of the stream as it approaches this point, resulting in an increase in the continuum optical depth, and hence a decrease of the Civ $`\lambda `$ 1550 Å equivalent width, and so is not necessarily indicative of a faint stream at this region.

### 6.2 Heating in the magnetically confined stream

Heating of the magnetic section of the stream has been modelled by Ferrario & Wehrse (1999) for a stream which is assumed to thread onto the field lines at a coupling radius $`r_c`$ (our $`R_\mu `$) from the white dwarf, over a radial distance $`\mathrm{\Delta }r_c`$ in the orbital plane. This provides an opportunity for a comparison between our observationally derived results, and the main points of their theoretical model results.

Ferrario & Wehrse (1999) consider two heating mechanisms: irradiation by the X-ray component from the accretion shock (the soft X-rays being up to $`4000`$ times more efficient at heating the stream than the hard X-rays) and the effects of magnetic reconnection in the stream-field interaction over the threading region $`\mathrm{\Delta }r_c`$. The comparison we make with these theoretical models is qualitative as the geometrical structure of the magnetically confined accretion flow in Ferrario & Wehrse is funnel-shaped, compared to our linear trajectory. However, our results showing brightness enhancements towards the threading region are consistent with their models which incorporated magnetic heating in the threading region. This implies that some magnetic heating mechanism is needed.

### 6.3 Temporal variations in the stream profile

Based on the colour ratios we deduced that the threading region was brighter and hotter in cycle 29995 than in cycle 29994 (Section 4.4). There are therefore significant changes in the threading region on the timescale of the orbital period (125 mins).

Although we have only modelled two cycles here, we can see that there is a difference in the brightness confirming the variability seen in the colour ratios and in the stream eclipse profiles from the light curves (see Figure 5). The large brightness enhancement near the white dwarf in cycle 29993 could be irradiated material which cools over the timescale of the next orbital period, leaving a fainter region as seen in cycle 29995 where the observed enhancement is less pronounced. However, without the model of the intermediate cycle, which we have excluded for reasons discussed previously, we cannot be certain and it will require more high signal-to-noise ratio observations of consecutive eclipses to investigate whether this in fact occurs. What is clear from the model results, and indeed directly from the raw light curves and colours, is that the emission from the whole stream is highly dynamic and unstable.

Dramatic changes in the stream eclipse profiles were also observed in the low state (Harrop-Allin et al. 2001, Glenn et al. 1994). These cycle-to-cycle changes in stream brightness and trajectory require a revision in our view of the stream, and in the manner in which the magnetic heating takes place in the threading region. Future treatments may need to consider large scale magnetic instabilities, and quasi-cyclic behavior.

## 7 Summary and conclusions

We have carried out high signal-to-noise ratio observations of HU Aqr using S-Cam2 on the WHT on two nights. The system was in a high accretion state, and from the single sharp change in the eclipse profile and archive soft X-ray light curves in this state, we infer that matter was accreting at only one pole on the white dwarf. At the onset of the eclipse, the accretion stream is the source of more than half of the optical emission from the system. The system brightness was similar from orbit to orbit, and the eclipse duration remained constant, but the shape of the accretion stream eclipse changed significantly, ending at $`\varphi =0.979`$, $`0.977`$ and $`0.982`$.

We find eclipse durations which are unchanged from past optical studies (Harrop-Allin et al. 1999b, Schwope et al. 1997), but shorter than those deduced in the soft X-rays by Schwope et al. (2001). However, the location of the accretion region is similar to that found by Schwope et al. (2001), so that there is no evidence that the optical emission is from higher latitudes on the white dwarf than the soft X-rays. This indicates that the duration of the soft X-ray eclipses may be affected by absorption, as they suggested.

The duration of the egress of the accretion region in the optical is $`8`$ s, compared to $`1.3`$ s in soft X-rays (Schwope et al. 2001). This indicates clearly that the region emitting cyclotron radiation is extended by a factor of $`5`$ by comparison with the soft X-ray emitting region, which Schwope et al. (2001) calculated as subtending an angle of $`3^{}`$.

We have found significant changes in the colour of the accretion stream from one eclipse to the next. This indicates that the threading region is hottest in the last of the eclipses by comparison with the previous two (see Section 4.4). We have modelled the stream using the technique of Hakala (1995) and Harrop-Allin et al. (1999a). This finds that most of the emission originates from two places, the region close to the white dwarf and in the threading region. By comparison with the models of Ferrario & Wehrse (1999) this indicates that magnetic heating is required in the threading region. The modelling clearly identifies an increase in brightness in the threading region for both cycles and an enhancement towards the white dwarf for cycle 29993. From this change in brightness of the heated regions in the model streams, and the varying stream eclipse profiles, we suggest that the magnetic heating in the threading region may be unstable. The implications of the highly variable stream trajectory and brightness profile should be recognised in future investigations of the stream properties. A systematic study of series of consecutive eclipses, with the highest possible signal-to-noise ratio, is required to investigate the characteristics of the magnetic heating and stream instabilities.

## ACKNOWLEDGMENTS

We acknowledge the contributions of other members of the Astrophysics Division of the European Space Agency at ESTEC involved in the optical STJ development effort, in particular S. Andersson, D. Martin, J. Page, P. Verhoeve, and J. Verveer. We acknowledge the excellent support given to the instrument’s operation at the WHT by the ING staff, in particular P. Moore and C.R. Benn.


text:
# Multiple Jet Production at Low Transverse Energies in 𝑝𝑝̄ Collisions at √𝑠 = 1.8 TeV

## I Introduction

The study of multiple jet production at high transverse energy was a goal of the 1993–1995 run of the Fermilab Tevatron collider, and the results have been compared with leading-order QCD predictions by both the CDF and DØ collaborations. These high-$`E_T`$ data, where $`E_T`$ is the transverse energy of the jet, are described satisfactorily by complete tree-level leading order $`2N`$ QCD calculations and by the herwig parton-shower Monte Carlo program. This kinematic region is described by $`Q^2/\widehat{s}1`$, where $`Q^2`$ is the square of the momentum transfer between partons (which we set equal to $`E_T^2`$), and $`\widehat{s}`$ is the square of the partonic center of mass energy. In this paper, we describe jet production measurements at significantly lower values of $`E_T`$ where detailed measurement of jet production in this kinematic region can provide information on the evolution of higher-order jet processes. In the same low $`E_T`$ region the DØ collaboration has previously reported the ratio of the inclusive three-jet to the inclusive two-jet cross section as a function of the scalar sum of jet transverse energies ($`H_T=E_T`$) with $`E_T>`$ 20 GeV. The ratio data can be described by the jetrad next-to-leading order Monte Carlo program. In this paper we make comparisons between Monte Carlo and several characteristics of multiple jet events including the leading jet transverse energy, the relative azimuthal angle between jets, and the summed vector transverse momenta of jets.

## II Data Sample and Corrections

The data were collected with the DØ detector at a proton-antiproton center-of-mass energy $`\sqrt{s}`$ = 1.8 TeV. Jets were identified using the liquid-argon uranium calorimeters, which have segmentation of $`\mathrm{\Delta }\eta \times \mathrm{\Delta }\varphi =0.1\times 0.1`$, where pseudorapidity $`\eta =\mathrm{ln}\mathrm{tan}\frac{\theta }{2}`$, $`\theta `$ is the polar angle, and $`\varphi `$ is azimuthal angle. At least one calorimeter trigger tower ($`\mathrm{\Delta }\eta \times \mathrm{\Delta }\varphi =0.2\times 0.2`$) with $`E_T`$ 2 GeV was required by a hardware trigger, and at least one jet with $`E_T`$ 12 GeV was required by a subsequent software trigger. Jets were reconstructed using a fixed cone algorithm with radius $`\mathrm{\Delta }=\sqrt{\mathrm{\Delta }\eta ^2+\mathrm{\Delta }\varphi ^2}=0.7`$ in $`\eta \varphi `$ space. The jet reconstruction threshold was $`E_T=8`$ GeV. If two jets overlapped and the shared transverse energy was more than 50% of the transverse energy of the lower-$`E_T`$ jet, the jets were merged; otherwise they were split into two jets. The integrated luminosity of this data sample is 2.0$`\pm `$0.3 nb-1. Instantaneous luminosity was restricted to be below 3 $`\times 10^{30}`$cm-2s-1 to minimize the number of multiple $`p\overline{p}`$ interactions in a single beam crossing.

To provide events of high quality, online and offline selection criteria suppressed multiple interactions, the cosmic ray background, and spurious jets. Jets were restricted to the pseudorapidity interval $`|\eta |3`$. The primary vertex of each event (reconstructed from time-of-flight as measured by scintillation counters) was required to be within 50 cm of the detector center.

Jet energies have been corrected for calorimeter response, shower development, and
various sources of noise. These corrections constitute the largest source of systematic uncertainty on the jet cross section. Typical values of the jet energy correction are (15–30)%, with an uncertainty of (2–4)%. In our study, we consider jets with $`E_T>`$ 20 GeV; for an inclusive $`n`$-jet event, the $`n`$ jets with the maximum $`E_T`$ (the leading jets) must have transverse energy above the threshold value. For example, a 3-jet event must have at least 3 jets above 20 GeV. The trigger efficiency is 0.85 for the inclusive ($`n`$ = 1) jet sample for energies near threshold, rising rapidly to unity at larger $`E_T`$. The efficiency is essentially unity for $`n>1`$.

To compare with data, Monte Carlo (MC) events were generated using the pythia 6.127, herwig 5.9, and jetrad programs. pythia and herwig simulate particle-level jets in the parton-shower approximation. jetrad simulates jets in the next-to-leading order approximation. To simulate detector resolution effects, the MC jet transverse energies were smeared with the experimentally determined jet energy resolution, which is $``$ 20% at $`E_T`$ = 20 GeV. Jet angular smearing used $`\eta `$ and $`\varphi `$ resolutions obtained by a MC simulation of the calorimeter response using herwig 5.9 and geant . These resolutions are $`0.08`$ at $`E_T`$ = 20 GeV. In pythia and herwig, jets were reconstructed at the particle level using the DØ algorithm, and in jetrad, at the parton level, using the Snowmass algorithm.

## III Leading jet $`E_T`$ distributions and systematic uncertainties

Distributions in transverse energy for the leading jet for inclusive $`n`$=1 to $`n`$=4 jet events are shown in Fig. 1, along with the results from pythia and herwig simulations. In these and all other plots, the data have been corrected for inefficiencies and energy calibration, but not for contributions from the underlying event. All simulated distributions have been smeared with energy and angular resolutions. Also to describe the data quantitatively, we normalize the theory (with a factor of 0.75 for pythia and 1.6 for herwig) to the observed two-jet inclusive cross section in Fig. 1(b) for $`E_T>40\mathrm{GeV}`$.

The normalised theory is in agreement with the data for all of the jet samples over the entire $`E_T`$ interval. A detailed comparison is shown in Figs. 2 and 3. Here the simulations have been brought into agreement with the data by selecting parameters that enhance low $`E_T`$ jet production. In the case of pythia, the core of the hadronic matter distribution has been increased to the fraction 0.32. An increased core fraction (the parameter PARP(83)) leads to enhancement of the multiple interaction rate, which tends to produce events with large multiplicity because of additional radiated low energy jets and underlying event energy. In the case of herwig, the minimum transverse momentum for the hard subprocesses has been set to 3.7 GeV. A decreased minimum transverse momentum (the parameter PTMIN) leads to increased soft underlying event contributions. The default values for these parameters are PARP(83)=0.5 and PTMIN=10$`\mathrm{GeV}`$. Variation of these values by more than 15% leads to disagreement with the low $`E_T`$ data. Other parameters, when varied from their default values, do not change the distributions significantly.

Figures 2 and 3 show the fractional difference (Data $``$ MC) / MC for the $`E_T`$ spectra in Fig. 1 with the uncertainties arising from jet-energy calibration and resolutions. The systematic uncertainty on the cross section is due primarily to the uncertainty in the energy calibration. This uncertainty can be estimated by considering cross sections derived with $`\pm 1`$ standard-deviation corrections to the jet energy scale. The same procedure can be used to derive the uncertainties due to jet energy and angular resolutions in the MC. At $`E_T`$ $`=`$ 25 GeV, the uncertainty in the three-jet cross section due to calibration of the data is 39%, and uncertainties in the MC due to energy and angular resolutions are 19% and 7%, respectively.

The uncertainty from energy resolution represents the dominant uncertainty in the MC. In Figs. 2 and 3, the relative systematic uncertainties in the cross section corresponding to the energy calibration added in quadrature with 15% uncertainty in luminosity are shown by the solid lines. The uncertainty in the ratio (Data $``$ MC) / MC from energy and angle smearing is shown by the dashed lines. The total uncertainty on the ratio is shown by the dotted lines. As shown in Figs. 2 and 3, both pythia and herwig describe the data quite well.

## IV Transverse energy and azimuthal distributions

To explore features of three- and four-jet production, we turn to observations of relative azimuthal distributions, distributions in summed transverse momenta, and three-jet studies. In Fig. 4(a) we plot the azimuthal difference between the leading two jets in events with two or more jets. Figures 4(b–d) show the azimuthal difference between the first and second, first and third, and second and third highest-$`E_T`$ jets in three-jet events. In Fig. 4(a) we see the strong anticorrelation (in the transverse plane) expected of two-jet events. The peak of the distribution widens substantially in the three-jet sample (Fig. 4(b–d)). The peaks correspond to the kinematic constraint of transverse momentum conservation for jets produced in hard QCD subprocesses. pythia (normalized as in Fig. 1) approximates the observed three-jet cross section and shapes. However, small discrepancies with herwig (also normalized as in Fig. 1) are evident.

Distributions of the square of the summed vector transverse momenta of jets $`Q_T^2=(𝐄_{T1}+𝐄_{T2}+\mathrm{}+𝐄_{Tn})^2`$ in Fig. 5 show significant imbalance of the transverse momenta for $`n`$ leading jets. If events at large $`Q_T^2`$ are removed by requiring balanced transverse energy, the corresponding three- and four-jet cross sections of Fig. 1 decrease at small $`E_T`$. The shoulder at $`Q_T^21600\mathrm{GeV}^2`$ in Fig. 5(a) can be eliminated by restricting the event sample to just two jets with $`E_T`$ above 20 GeV, and no other jets between 8 and 20 GeV. This shoulder can consequently be associated with higher-order radiation.

To find the pair of jets $`\{i,j\}`$ most likely to originate from the hard interaction (rather than from gluon bremsshtrahlung), we define the scaled summed dijet vector transverse momentum: $`𝐪_{ij}=(𝐄_{Ti}+𝐄_{Tj})/(E_{Ti}+E_{Tj})`$. We choose the pair with the smallest magnitude of this vector and in Figs. 6(a) and 7(a) plot the distribution of the relative azimuthal angle $`\mathrm{\Phi }_c`$ between the jets in that pair. The data, pythia, and herwig show a narrow maximum in the region where two jets from the hard scatter appear back-to-back ($`\mathrm{\Phi }_c=\pi `$). The prediction from jetrad is peaked away from $`\mathrm{\Phi }_c\pi `$ because only one extra jet is present.

Figures 6(b,c) and Fig. 7(b,c) show the azimuthal separation of the third jet from each of the two jets that correspond to the minimum $`q_{ij}^2`$. These distributions contain events only for $`\pi \mathrm{\Phi }_c0.4`$; that is, events in which the balanced jets are essentially back-to-back. If the third jet were correlated with the balanced jets, it would be observed nearby or opposite the balanced jets. However, the data show the third jet to be weakly correlated with the balanced jets, and emitted at all angles. The uncertainties associated with energy calibration and luminosity are shown by the solid lines in Figs. 6 and 7. Uncertainties from the energy resolution are shown by dashed lines in Fig. 6.

We see that the data, pythia, and herwig have wider distributions than jetrad. pythia describes the data quite well, while jetrad fails. The agreement with pythia has been achieved only with enhanced multiple parton interaction rates. herwig demonstrates small qualitative disagreement with the shape of the azimuthal plot of Fig. 7(b); the peak at $`\pi /2`$ is produced by jets reconstructed from the underlying-event energy and grows quickly with small changes in PTMIN. Such jets are strongly overlapped with more that one jet. If jets overlapping two or more nearby jets are excluded, the herwig shape in Fig. 7(b) improves but the agreement shown in Fig. 7(a) worsens. (The cone algorithm reconstructs jets from seed towers and may therefore reconstruct jets sharing energy. The reconstruction algorithm then merges or splits the energy encompassed in these overlapping jets.) Elimination of these jets tends to suppress contributions from the soft underlying event. Soft interactions result in a wide distribution of particles throughout angular phase space. Jets reconstructed from these particles tend to be wider and of lower energy than more collimated partonic jets. Such jets often share a significant fraction of energy with similar, neighbouring jets and are merged into a single jet.

The shapes of the simulated distributions are sensitive to modeling of the multiple parton interactions. Tuning of the multiple interaction contribution in pythia and the minimum generated transverse momentum in herwig are required for good agreement. In particular, simulations with smaller contributions from soft parton interactions show discrepancies with the data.

## V Conclusions

In this paper we showed comparisons between Monte Carlo and data for several characteristics of multiple jet events with a low jet-$`E_T`$ threshold. These comparisons included the leading jet transverse energy, the relative azimuthal angle between jets, and the summed vector transverse momenta of jets. Our data on multiple jet production at low $`E_T`$ agree with pythia and herwig. This is observed in the distributions of the transverse energy of the leading jets (Fig. 1), azimuthal distributions (Fig. 4), in the square of the summed vector transverse momenta $`Q_T^2`$ (Fig. 5), and in the three-jet angular distributions that suggest the presence of a weakly correlated jet (Figs. 6, 7). jetrad cannot adequately describe the angular distributions of the three leading jets in three jet events.

## VI Acknowledgments

We thank the staffs at Fermilab and collaborating institutions, and acknowledge support from the Department of Energy and National Science Foundation (USA), Commissariat à L’Energie Atomique and CNRS/Institut National de Physique Nucléaire et de Physique des Particules (France), Ministry for Science and Technology and Ministry for Atomic Energy (Russia), CAPES and CNPq (Brazil), Departments of Atomic Energy and Science and Education (India), Colciencias (Colombia), CONACyT (Mexico), Ministry of Education and KOSEF (Korea), CONICET and UBACyT (Argentina), The Foundation for Fundamental Research on Matter (The Netherlands), PPARC (United Kingdom), Ministry of Education (Czech Republic), A.P. Sloan Foundation, and the Research Corporation.


text:
# SCALE-FREE EQUILIBRIA OF SELF-GRAVITATING GASEOUS DISKS WITH FLAT ROTATION CURVES

## 1 INTRODUCTION

Astrophysical disks are either stellar or gaseous. The existence of non-axisymmetric equilibrium states of such systems has been attractive both from theoretical and observational aspects. The search for dynamical equilibria of stellar systems is known as the self-consistency problem where one seeks for a positive phase space distribution function (DF) that produces the density of the model and satisfies the collisionless Boltzmann equation (Binney & Tremaine 1987). The construction of the DF can be done by analytical means or numerical exploration based on Schwarzschild’s orbit superposition method (Schwarzschild 1979; Kuijken 1993, hereafter K93; Levine & Sparke 1998; Teuben 1987; de Zeeuw, Hunter, & Schwarzschild 1987; Jalali & de Zeeuw 2002, hereafter JZ). The main difficulty linked to the self-consistency problem is the lack of the second integral of motion for most planar models.

Gaseous disks, however, are investigated through solving the hydrodynamic equations of compressible, self-gravitating fluids. There are several works in the literature that deal with the stability of perturbed axisymmetric gaseous disks (Heemskerk, Papaloizou, & Savonije 1992; Goodman & Evans 1999; Lemos, Kalnajs, & Lynden-Bell 1991). The problem of non-axisymmetric systems is somewhat different because our knowledge is limited even about equilibrium states. Following a first-order linear theory, Syer & Tremaine (1996, hereafter ST96) constructed simple scale-free disks and approximately solved the steady-state hydrodynamic equations for the components of the velocity vector. They also proposed a nonlinear theory based on the calculation of Fourier coefficients using Newton’s method. A similar strategy was adopted by Galli et al. (2001, hereafter G01) in the study of isothermal molecular clouds and the fragmentation hypothesis of multiple stars. Their nonlinear solutions, in contrast with ST96, were generated by direct numerical integration of reduced ordinary differential equations (ODE) and an iterative scheme to match the boundary conditions. These methods of approach seem inefficient when one seeks for global properties of highly non-axisymmetric equilibria and the transient response of fluid disks in the vicinity of such states. Therefore, finding analytical solutions for the hydrodynamic equations is of special interest even in the steady-state conditions.

In this paper we explore exact non-axisymmetric solutions of self-similar gaseous disks with flat rotation curves. In §2 we express the governing equations and reduce them to a third-order ODE in §3 where we also introduce our closed-form analytic solutions. We construct hydrodynamical counterparts of the JZ models in §4 and apply our results to the Large Magellanic Cloud (LMC) in §5. We show that the kinematics of the LMC can be explained by our model elliptic disks. We end up with conclusions in §6.

## 2 STEADY-STATE HYDRODYNAMIC EQUATIONS

We consider a thin disk of compressible, inviscid self-gravitating fluid in a steady state (dynamical equilibrium). The hydrodynamic equations can be written as the time-independent continuity equation

$$\frac{1}{R}\frac{}{R}(R\mathrm{\Sigma }v_R)+\frac{1}{R}\frac{}{\varphi }(\mathrm{\Sigma }v_\varphi )=0,$$
(1)
and the momentum equations

$`{\displaystyle \frac{1}{R}}{\displaystyle \frac{}{R}}(R\mathrm{\Sigma }v_R^2)+{\displaystyle \frac{1}{R}}{\displaystyle \frac{}{\varphi }}(\mathrm{\Sigma }v_Rv_\varphi ){\displaystyle \frac{\mathrm{\Sigma }v_\varphi ^2}{R}}={\displaystyle \frac{p}{R}}\mathrm{\Sigma }{\displaystyle \frac{U}{R}},`$ (2)
$`{\displaystyle \frac{1}{R}}{\displaystyle \frac{}{R}}(R\mathrm{\Sigma }v_Rv_\varphi )+{\displaystyle \frac{1}{R}}{\displaystyle \frac{}{\varphi }}(\mathrm{\Sigma }v_\varphi ^2)+{\displaystyle \frac{\mathrm{\Sigma }v_Rv_\varphi }{R}}={\displaystyle \frac{1}{R}}{\displaystyle \frac{p}{\varphi }}\mathrm{\Sigma }{\displaystyle \frac{1}{R}}{\displaystyle \frac{U}{\varphi }},`$ (3)
where $`(R,\varphi )`$ are the polar coordinates, $`𝐯=(v_R,v_\varphi )`$ is the velocity field in Eulerian description, and $`p`$ is the pressure. The potential $`U`$ is due to the self-gravity, which is related to the surface density $`\mathrm{\Sigma }`$ through the Poisson integral

$$U(R,\varphi )=G\frac{\mathrm{\Sigma }(R^{},\varphi ^{})R^{}\mathrm{d}R^{}\mathrm{d}\varphi ^{}}{\sqrt{R^2+R^22RR^{}\mathrm{cos}(\varphi \varphi ^{})}},$$
(4)
where $`G`$ is the gravitational constant. In the following section, we attempt to simplify the governing equations for scale-free disks and solve the resulting equations by analytical means.

## 3 NON-AXISYMMETRIC DISKS WITH LOGARITHMIC POTENTIALS

Gaseous disks with logarithmic potentials possess the astrophysically important property of a flat rotation curve. The circular gaseous disks having this property are known as Mestel (1963) disks, the stability of which was investigated by Goodman & Evans (1999). Here we derive exact analytical expressions for the physical quantities of non-axisymmetric disks.

### 3.1 Stream Function

We assume cuspy, scale-free potential–density pairs of the form

$`U(R,\varphi )`$ $`=`$ $`U_0[2\mathrm{ln}R+g(\varphi )],`$
$`\mathrm{\Sigma }(R,\varphi )`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Sigma }_0f(\varphi )}{R}},`$ (5)
where $`U_0`$ and $`\mathrm{\Sigma }_0`$ are constants. The functions $`f(\varphi )`$ and $`g(\varphi )`$ are $`2\pi `$-periodic in $`\varphi `$ because they are continuous. We normalize the surface density so that

$$\frac{1}{2\pi }_\pi ^\pi f(\varphi )d\varphi =1.$$
(6)
This transfers the scaling of $`\mathrm{\Sigma }`$ to $`\mathrm{\Sigma }_0`$ (e.g., G01) and implies $`U_0=\pi G\mathrm{\Sigma }_0`$. Given $`f`$, Poisson’s integral does not usually lead to a closed-form expression for $`g`$ and one has to adopt a series solution. Let us assume that $`f(\varphi )`$ admits the Fourier expansion

$$f(\varphi )=1+\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\left(C_{\mathrm{}}\mathrm{cos}\mathrm{}\varphi +S_{\mathrm{}}\mathrm{sin}\mathrm{}\varphi \right),$$
(7)
where $`C_{\mathrm{}}`$ and $`S_{\mathrm{}}`$ are constant coefficients. Using eqs (A1) and (A5) of ST96, we find the potentials corresponding to single Fourier modes and superpose the results. This gives

$$g(\varphi )=2\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\left(\frac{C_{\mathrm{}}}{\mathrm{}}\mathrm{cos}\mathrm{}\varphi +\frac{S_{\mathrm{}}}{\mathrm{}}\mathrm{sin}\mathrm{}\varphi \right).$$
(8)
We prescribe $`f`$ and determine $`g`$ through the scheme mentioned above. Then, we attempt to solve the momentum equations for $`v_R`$, $`v_\varphi `$, and the pressure $`p`$ in terms of $`f`$ and $`g`$.

The velocity components can be derived from the stream function $`\mathrm{\Psi }`$ as

$$\mathrm{\Sigma }v_R=\frac{1}{R}\frac{\mathrm{\Psi }}{\varphi },\mathrm{\Sigma }v_\varphi =\frac{\mathrm{\Psi }}{R}.$$
(9)
The scale-freeness of the model suggests the Ansatz

$$\mathrm{\Psi }(R,\varphi )=\mathrm{\Sigma }_0V_\varphi \left[\mathrm{ln}R+P(\varphi )\right],$$
(10)
where $`V_\varphi `$ is a constant reference velocity that will become the circular velocity if the disk is axisymmetric ($`f=1`$, $`P=0`$). We utilize eqs (5) through (10) and take the curl of eqs (2) and (3). This eliminates the pressure $`p`$ and leaves us with the following ODE for $`P(\varphi )`$

$$P^{\prime \prime \prime }\frac{2f^{}}{f}P^{\prime \prime }+\frac{1}{f^2}\left[2(f^{})^2ff^{\prime \prime }+f^2\right]P^{}=\beta f\left(fg^{}+2f^{}\right),$$
(11)
where $`\beta =U_0/V_\varphi ^2`$ is a dimensionless parameter and the prime denotes $`\mathrm{d}/\mathrm{d}\varphi `$. The solution $`P(\varphi )`$ will result in the velocity field of our planar self-gravitating fluids. Physical solutions of (11) correspond to the condition $`𝐯(R,\varphi )=𝐯(R,\varphi +2\pi )`$, which ensures the continuity of velocity field. Hence, we seek for $`2\pi `$-periodic solutions of $`P(\varphi )`$.

### 3.2 Velocity Components

Carrying out a change of dependent variable as (this is indeed Liouville’s transformation explained in Appendix A)

$$P^{}=w(\varphi )\frac{f(\varphi )}{f(0)},$$
(12)
eq. (11) reads

$$w^{\prime \prime }+w=\beta f(0)F_0(\varphi ),F_0=fg^{}+2f^{}.$$
(13)
The particular solution of (13) is then obtained as

$$w(\varphi )=\beta f(0)\left[a_0\mathrm{sin}\varphi +b_0\mathrm{cos}\varphi +_0^\varphi F_0(s)\mathrm{sin}(\varphi s)ds\right],$$
(14)
where the (dimensionless) constant coefficients $`a_0`$ and $`b_0`$ are given by

$$\{\begin{array}{c}a_0\\ b_0\end{array}\}=\frac{1}{\pi }_0^{2\pi }\{\begin{array}{c}\mathrm{sin}\varphi \\ \mathrm{cos}\varphi \end{array}\}d\varphi _0^\varphi F_0(s)\mathrm{sin}(\varphi s)ds.$$
(15)
The function $`w(\varphi )`$, and therefore $`𝐯`$, will be $`2\pi `$-periodic if we impose the boundary conditions $`w(\varphi +2\pi )=w(\varphi )`$ and $`w^{}(\varphi +2\pi )=w^{}(\varphi )`$ that lead to the necessary conditions

$$_0^{2\pi }F_0(\varphi )\mathrm{sin}\varphi \mathrm{d}\varphi =_0^{2\pi }F_0(\varphi )\mathrm{cos}\varphi \mathrm{d}\varphi =0.$$
(16)
Using (12) and (9), the velocity components are readily determined as

$`v_R`$ $`=`$ $`\beta V_\varphi \left[a_0\mathrm{sin}\varphi +b_0\mathrm{cos}\varphi +{\displaystyle _0^\varphi }F_0(s)\mathrm{sin}(\varphi s)ds\right],`$
$`v_\varphi `$ $`=`$ $`{\displaystyle \frac{V_\varphi }{f(\varphi )}},`$ (17)
from which one can obtain $`P(\varphi )`$. It is now straightforward to determine the pressure $`p`$ from either (2) or (3) as

$$p=\frac{\mathrm{\Sigma }_0V_\varphi ^2}{Rf(\varphi )}\left[\beta A(\varphi )1\right]=c_s^2\mathrm{\Sigma },$$
(18)
where

$$c_s^2=\frac{V_\varphi ^2}{f(\varphi )^2}[\beta A(\varphi )1],$$
(19)
is the square of the sound speed (at the azimuth $`\varphi `$) and

$$A(\varphi )=f(\varphi )\left[(b_0\mathrm{sin}\varphi a_0\mathrm{cos}\varphi )+2f(\varphi )_0^\varphi F_0(s)\mathrm{cos}(\varphi s)ds\right].$$
(20)
Sound waves are stable if $`c_s^2>0`$. i.e., $`p>0`$. Define $`A_{\mathrm{min}}=\mathrm{min}\left\{A(\varphi ):0\varphi 2\pi \right\}`$. The necessary and sufficient condition for the positiveness of $`p`$ will be

$$0<\frac{1}{\beta }<A_{\mathrm{min}}.$$
(21)
Apart from constraining the admissible values of $`\beta `$, this condition may put some restrictions on the model parameters. As it can be seen in (3.2), $`v_\varphi `$ is maximum (minimum) where the surface density is minimum (maximum). This result has previously been used in the construction of DFs of stellar systems using Schwarzschild’s orbit superposition method (e.g., K93; JZ).

Our scale-free analytical solutions are general. They are independent of special angular variation of surface density and its corresponding potential. Our solutions can be used even in the presence of an external force that falls off as steep as $`R^1`$. Such a force can be that of a magnetic field or the gravitational influence of a dark matter halo. These effects can enter in our equations through $`F_0(\varphi )`$.

The most important point related to our solutions is the existence of transonic flows. According to (19), it is possible to obtain entirely subsonic, transonic, and entirely supersonic flows depending on the choice of $`\eta =1/\sqrt{\beta A_{\mathrm{min}}}`$ in the range $`0<\eta <1`$. Let us define the Mach number

$$=\frac{1}{c_s}\sqrt{v_R^2+v_\varphi ^2},$$
(22)
which is a function of the azimuthal angle $`\varphi `$ in our disks. For convenience, we use the variable

$$_0(\varphi )=\frac{1}{+1},$$
(23)
for the classification of flows. Disks with $`_0(\varphi )>0`$ ($`_0(\varphi )<0`$) for $`\varphi [0,2\pi ]`$ are entirely supersonic (subsonic) and our steady-state models are credible. Flows for which $`_0(\varphi )`$ switches sign along a streamline, are transonic. In such a circumstance, steady-state models become problematic as explained below. In principle, transition from subsonic to supersonic speeds (along a streamline) can take place smoothly, but the reverse is impossible without assuming a shock wave. Because fluid elements in the subsonic part of a streamline, have not information from the arrival of the particles in the supersonic upstream. Therefore, they suddenly find themselves collided with these particles. After a while, these interactions develop a shock wave and the speed of supersonic stream discontinuously decreases to a subsonic value. This process is inherently time-dependent and it increases entropy and generates heat. Thus, we must cautiously deal with any transonic solution. Galli et al. (2001) also argue that models with transonic flows will most likely be unstable. Despite these arguments, there is an exceptional case in nature where a supersonic flowfield is isentropically (and of course smoothly) compressed from supersonic to subsonic speeds. That is the famous conical flow detailed description of which can be found in Anderson (1990, chapter 10). In this very special case, again, a shock wave still presents but the sonic line does not coincide with the shock ray. The latter completely lies in the supersonic region.

In any case, it seems reasonable to exclude transonic solutions from our results. A stability analysis will surely provide much insight into the problem. The instabilities of our non-axisymmetric disks will be resolved in a future paper.

For axisymmetric disks, the parameter $`\eta `$ is obtained in terms of Mach number through

$$\eta =\frac{}{\sqrt{1+^2}}.$$
(24)
In this case, $``$ is independent of $`\varphi `$ and the transition from subsonic to supersonic flows occurs at $`\eta =1/\sqrt{2}`$.

## 4 EXAMPLES

We now apply our analytical tool to certain non-axisymmetric models for which we have either observational evidence or stellar dynamical counterparts.

### 4.1 Simple Disks

The simplest form of non-axisymmetric disks is obtained when a circular disk is perturbed by a single Fourier mode as

$$f(\varphi )=1+f_m\mathrm{cos}m\varphi ,$$
(25)
where $`0|f_m|<1`$. The corresponding $`\varphi `$-dependence of the potential will become

$$g(\varphi )=\frac{2}{m}f_m\mathrm{cos}m\varphi .$$
(26)
One can readily verify that the condition (16) is satisfied for all $`m1`$ and from (3.2) we find

$$v_R=\beta V_\varphi f_m\left(\frac{f_m}{4m^21}\mathrm{sin}2m\varphi \frac{2}{m+1}\mathrm{sin}m\varphi \right).$$
(27)
It is not difficult to show that $`A(\varphi )`$ is positive for our simple disks when $`0|f_m|<1`$. These disks are analogous to the multiple-lobed disks of G01. But they are not the same. The difference can be traced back to the logic of solving the hydrodynamic equations. In fact, hydrodynamic equations for self-gravitating disks are four equations for the five unknowns $`\mathrm{\Sigma }`$, $`U`$, $`v_R`$, $`v_\varphi `$ and $`p`$. Thus, we are free to specify one of these variables, or put an extra constraint on them. Galli et al. (2001) follow the latter approach and assume an equation of state in the form $`p=a^2\mathrm{\Sigma }`$ where $`a`$ is the constant sound speed. They numerically solve their problem for $`\mathrm{\Sigma }`$, $`v_\varphi `$, and $`v_R`$, and obtain $`U`$ through Poisson’s equation. Our disks, however, are constructed relying on another approach. We prescribe the surface density $`\mathrm{\Sigma }`$ and find $`U`$ using Poisson’s integral. Then, we determine $`v_R`$, $`v_\varphi `$, and $`p`$. Our disks, in contrast to G01, are not isothermal \[because the sound speed is not constant according to (19)\]. Nonetheless, our solutions are exact and we never face the convergence problem as G01.

In Figure 1, we have surveyed the $`(\eta ,f_m)`$ parameter space of our simple disks for $`\pi G=1`$, $`\mathrm{\Sigma }_0=1`$, and several choices of the wavenumber $`m`$. For each $`m`$, we obtain a curve (showed by solid, dashed, and dash-dotted lines) which is the boundary of models with transonic flowfields. The region to the up side of each curve (indicated by the letter B) corresponds to the models that are transonic along a typical streamline. Below the boundary curve, we obtain entirely subsonic and entirely supersonic regions. Excluding transonic solutions (they are unphysical due to arguments in §3.2), highly non-axisymmetric disks can be constructed in the subsonic region as the wavenumber $`m`$ is increased. The maximum elongation of the lopsided mode ($`m=1`$) occurs for $`f_m0.5`$ and $`\eta 1/\sqrt{2}`$. Whatever the value of $`m`$ may be, the maximum value of $`\eta `$ is $`1/\sqrt{2}`$ for entirely subsonic models. Entirely supersonic models are approximately circular with $`f_m0.2`$.

### 4.2 Bisymmetric Disks

There are two interesting families of bisymmetric disks for which self-consistent stellar dynamical models are known. They are the power-law and elliptic disks (JZ). Disks whose equipotentials are similar concentric ellipses, are called the power-law disks. They are planar analogs of the logarithmic ellipsoidal models (e.g., Binney 1981; Evans, Carollo & de Zeeuw 2000). The power-law disks with flat rotation curves take the potential of the form (5) with

$$g(\varphi )=\mathrm{ln}(q^2\mathrm{cos}^2\varphi +\mathrm{sin}^2\varphi ),$$
(28)
where $`q`$ is the axial ratio of equipotentials. The angular dependence of the surface density corresponding to (28) has been derived in JZ as the following convergent series

$`f(\varphi )`$ $`=`$ $`1+{\displaystyle \underset{\mathrm{even}\mathrm{}>0}{\overset{\mathrm{}}{}}}f_{\mathrm{}}\mathrm{cos}\mathrm{}\varphi `$ (29)
$`f_{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left\{{\displaystyle \frac{(1q)^n\delta _k\mathrm{}}{2^{n1}}}+{\displaystyle \underset{m=0}{\overset{n1}{}}}{\displaystyle \frac{(1q)^n\beta _m^{(n)}}{2n}}\left[i\delta _i\mathrm{}+j\delta _j\mathrm{}\right]\right\},`$
with $`k=2n`$, $`i=n+m`$, and $`j=nm`$. In this relation, $`\delta _{rs}`$ is the Kronecker delta and we have defined $`\beta _m^{(n)}=0`$ if $`nm`$ is odd, and

$`\beta _m^{(n)}`$ $`=`$ $`\{\begin{array}{cc}\frac{1}{2^{n1}}(\begin{array}{c}n\\ {\displaystyle \frac{nm}{2}}\end{array}),\hfill & m0,\hfill \\ & \\ \frac{1}{2^n}(\begin{array}{c}n\\ n/2\end{array}),\hfill & m=0,\hfill \end{array}`$ (37)
if $`nm`$ is even.

When the isodensity contours are similar concentric ellipses, we attain elliptic disks with

$$f(\varphi )=(q^2\mathrm{cos}^2\varphi +\mathrm{sin}^2\varphi )^{1/2},$$
(38)
where $`q`$ is the axis ratio of isodensity contours. Stellar elliptic disks were constructed (for the first time) in K93 and then were examined for self-consistency by the curvature condition in JZ. The potential function associated with (38) was determined by Evans & de Zeeuw (1992, eq. 5.2) in terms of an indefinite integral, but a practical way is to generate it using the scheme presented in §3.1. We find

$`U(R,\varphi )`$ $`=`$ $`2\pi G\mathrm{\Sigma }_0\left(\mathrm{ln}R{\displaystyle \underset{\mathrm{even}\mathrm{}>0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{f_{\mathrm{}}}{\mathrm{}}}\mathrm{cos}\mathrm{}\varphi \right),`$ (39)
$`f_{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \underset{m=0}{\overset{n}{}}}{\displaystyle \frac{(1q)^n\mathrm{\Gamma }(\frac{1}{2}+k)\mathrm{\Gamma }(\frac{1}{2}+nk)\beta _m^{(n)}(\delta _i\mathrm{}+\delta _j\mathrm{})}{2\pi k!(nk)!}},`$
with $`i=2kn+m`$ and $`j=2knm`$. All the series introduced so far, rapidly converge using only the first few terms.

For instance, we have computed the isodensity contours and streamlines of our bisymmetric disks for $`q=0.5`$, $`\pi G=1`$, and $`\mathrm{\Sigma }_0=1`$, which result in $`A_{\mathrm{min}}=0.854`$ and $`A_{\mathrm{min}}=3.041`$ for the power-law and elliptic disks, respectively. Figure 2 shows the results of our calculations for several choices of $`V_\varphi `$. Solid lines are the isodensity contours and dashed lines are streamlines. For the small values of $`V_\varphi `$, the pressure increases and dominates centrifugal forces. Therefore, streamlines become highly dimpled in both disks (Panels 2a and 2c). As the value of $`V_\varphi `$ is increased ($`V_\varphi ^2U_0A_{\mathrm{min}}`$, i.e., $`\beta 1/A_{\mathrm{min}}`$), centrifugal forces overwhelm the pressure forces and keep the fluid particles away from the center. This causes streamlines to become rounder (Panels 2b and 2d). Numerical computations show that $`A(\varphi )`$ ($`0\varphi 2\pi `$) is positive for the power-law and elliptic disks in the range $`0q1`$. Therefore, it is always possible to find a physical value for $`V_\varphi `$. Figure 3 shows the $`(\eta ,q)`$ parameter space of the power-law and elliptic disks. The region filled by dots corresponds to models with transonic flowfields. They are most likely unphysical models as we discussed in §3.2. Hence, we exclude them from our solution space. As Figure 3 shows, supersonic elliptic disks can be more elongated than supersonic power-law disks. In both cases the maximum value of $`\eta `$ for subsonic flows is equal to $`1/\sqrt{2}`$. This is a generic property of our self-gravitating gaseous disks inherited from simple disks.

Another remarkable point is the limit of the axis ratio $`q`$ for plausible models (transonic models excluded). The lower bound of $`q`$ is almost 0.4 in elliptic disks but in the case of power-law disks one can see a narrow subsonic region with $`q<0.08`$. In the limit of $`q0`$, these highly elongated models converge to a needle. The results for elliptic disks are similar to the findings of JZ for stellar disks. In fact, stellar elliptic disks with $`q0.5`$ are non-self-consistent. On the evidence of Figure 3a, gaseous power-law disks exhibit significantly different features compared to their stellar counterparts that cannot be self-consistent if $`q<0.707`$ (JZ).

Our calculations show that the gas pressure of elliptic disks approximately takes an axisymmetric distribution when $`V_\varphi `$ is small. As a result, the line-of-sight velocity dispersion $`\sigma =v_{\mathrm{los}}^2v_{\mathrm{los}}^2`$ will inversely be proportional to $`f(\varphi )`$. This property does not depend on the choice of $`q`$ and it predicts a correlation between the photometric and spectroscopic data of elliptic disks when centrifugal forces are dominated by internal stresses.

## 5 IMPLICATIONS FOR THE KINEMATICS OF THE LMC

In his recent study of the LMC, van der Marel (2001) has shown that the LMC disk is elongated and its isophotal lines are well fitted by ellipses. Since the LMC has an approximately flat rotation curve, our elliptic disks (presented in §4.2) can be applied to model the LMC disk.

The most recent studies of the kinematics of gas and discrete tracers (like carbon stars) of the LMC (Kim et al. 1998; Alves & Nelson 2000) show that the line of maximum velocity gradient is in the range $`143^{}<\mathrm{\Theta }_{\mathrm{max}}<183^{}`$, which differs from the angular position of the line of nodes ($`\mathrm{\Theta }=122^{}.5`$) by $`20^{}<\mathrm{\Theta }_{\mathrm{max}}\mathrm{\Theta }<60^{}`$. Noting the fact that the intrinsic ellipticity of the LMC disk is $`ϵ=0.312`$ in outer regions, van der Marel (2001) argues that $`\mathrm{\Delta }\mathrm{\Theta }=\mathrm{\Theta }_{\mathrm{max}}\mathrm{\Theta }20^{}`$ is possible but it is hard to believe that $`\mathrm{\Delta }\mathrm{\Theta }`$ can be as large as $`60^{}`$. He ascribes the excess of computed misalignment to the center-of-mass and rigid-body motions of the LMC. A more accurate kinematic analysis of the LMC predicts a value of $`\mathrm{\Delta }\mathrm{\Theta }`$ as small as $`7.4^{}`$ (van der Marel 2002, private communication). To provide a theoretical framework for these observations, we have modeled the outer region of the LMC using our elliptic disks.

We have calculated the position angle of the maximum line-of-sight velocity of our elliptic disks, given the viewing angles $`\mathrm{\Theta }=122^{}.5`$ and $`i=34^{}.7`$ that were determined by van der Marel & Cioni (2001). Our results have been shown in Figure 2. This figure illustrates the variation of $`\mathrm{\Delta }\mathrm{\Theta }`$ versus $`\eta `$ for several choices of $`ϵ=1q`$ including $`ϵ=0.312`$, which is the intrinsic ellipticity of the LMC at large radii. A discontinuity occurs in the graph of $`\mathrm{\Delta }\mathrm{\Theta }`$ (let us say at $`\eta =\eta _{\mathrm{cr}}`$) because of a sudden change in the morphology of streamlines. For $`\eta <\eta _{\mathrm{cr}}`$, streamlines are dimpled. The graph for $`ϵ=0.312`$ (filled squares) shows that $`\mathrm{\Delta }\mathrm{\Theta }`$ varies between $`0^{}`$ and $`30^{}`$ when streamlines are not dimpled. Interestingly, this result is in harmony with the predictions of van der Marel (2001). It is remarked that for $`ϵ=0.312`$ we obtain $`\eta _{\mathrm{cr}}=0.36`$. As $`\eta 1`$, streamlines become almost circular independent of the value of $`ϵ`$. In such a circumstance, the position angle of the maximum line-of-sight velocity coincides with that of the line of nodes, i.e. $`\mathrm{\Delta }\mathrm{\Theta }0`$, as expected. Based on the available observational data and our results for $`ϵ=0.312`$, we conclude that our scale-free elliptic disks (with flat rotation curves) can explain the kinematics of the outer region of the LMC if $`\eta `$ is chosen in the range $`0.36<\eta <0.72`$. The upper bound on $`\eta `$ comes from the fact that our transonic solutions are not reliable answers of steady-state equations. For $`\mathrm{\Delta }\mathrm{\Theta }=7.4^{}`$, we find $`\eta 0.54`$, which is well within the acceptable range. Moreover, models with $`\eta 1`$ (they are supersonic) lead to very high fluctuations of the pressure, and therefore, velocity dispersion. This does not seem to be relevant to a stable model. We also rule out the possibility of $`\eta <\eta _{\mathrm{cr}}`$, for it gives $`\mathrm{\Delta }\mathrm{\Theta }<0`$, which is not observed. Furthermore, the condition $`p0`$ at the boundary of the LMC disk is inconsistent with a pressure-dominated model.

## 6 CONCLUSIONS

Exact equilibrium solutions of self-gravitating systems are rare. Two classes of such solutions were found by Medvedev & Narayan (2000) for three-dimensional, axisymmetric isothermal systems. A generalization of these solutions for non-axisymmetric rotating systems was then introduced by Shadmehri & Ghanbari (2001). Although the density distribution of their models is non-axisymmetric, the velocity field has only the toroidal component. This makes their models too simple to be used in the study of fully three-dimensional non-axisymmetric distributions of matter.

In a systematic search for non-axisymmetric equilibrium solutions, we started our study from gaseous disks because they inherit many properties of three-dimensional systems. Using the benefit of the existence of a stream function, we were successful in reducing the governing equations of scale-free models to an incomplete linear ODE for the angular part of the stream function. It then became clear that Liouville’s transformation can reduce the obtained ODE to a second-order inhomogeneous ODE, which is similar to that of a forced harmonic oscillator. This was striking because it led us to derive closed-form solutions for the radial velocity component, stream function and the pressure $`p`$. In fluid mechanics, transport equations are written for mass, momentum and energy. While the pressure occurs in governing equations, there is no extra independent equation in order to determine $`p`$. The problem is rather difficult in the presence of an anisotropic stress tensor, which is encountered in turbulence and Jeans’ equations of stellar hydrodynamics. Taking the curl of momentum equations is a trick that temporarily removes the pressure from our calculations of the stream function. But, one should take care of the sign of $`p`$ once it is eventually extracted from momentum equations. Only a positive pressure distribution is physical because a fluid medium can resist against a compressive stress not a tensional one.

Non-axisymmetric stellar disks have necessarily anisotropic DFs and pressure tensors. This point is the main difference between the physics of a steady-state gaseous disk and its stellar counterpart described by Jeans’ equations. Nevertheless, isotropic pressure distribution is not an obstacle for the existence of non-axisymmetric, self-gravitating gaseous disks as we showed in §4.

Our scale-free elliptic disks were used to model the outer part of the LMC disk. The main restriction of our modeling is the isotropy of the stress tensor. We speculate that our results will not be changed in a significant way if the model is stellar with an anisotropic pressure. The other important issue is the perturbation of the Galaxy. The tidal field of the Galaxy induces a force whose $`R`$-distribution is totally different from that of a logarithmic potential field. However, if the tidal field of the Galaxy was dominant, the rotation curve of the LMC would show substantial deviation from flatness. But, we know that the rotation curve of the LMC is almost flat. Therefore, our isolated models provide acceptable zeroth-order solutions. Such a modeling of the LMC is a quick route for interpreting the observational data. It is obvious that taking the tidal effects into account, will destroy the enjoyable scale-free nature of our disks.

We thank Tim de Zeeuw for useful comments and drawing our attention to the possible application of our findings to the LMC. We are indebted to an anonymous referee for enlightening suggestions that helped us to substantially improve the presentation of our results.

## Appendix A LIOUVILLE’S TRANSFORMATION

Consider the following non-autonomous, linear, ordinary differential equation for $`Z`$ in terms of $`t`$:

$$\frac{\mathrm{d}^2Z}{\mathrm{d}t^2}+Q_1(t)\frac{\mathrm{d}Z}{\mathrm{d}t}+Q_2(t)Z=Q_0(t).$$
(A1)
Using the transformation (Bellman 1966, page 82; Verhulst 1996, problem 8-2)

$$Z(t)=w(t)\mathrm{exp}\left(\frac{1}{2}_0^tQ_1(\tau )d\tau \right),$$
(A2)
which has been introduced by Liouville, eq. (A1) becomes

$$\frac{\mathrm{d}^2w}{\mathrm{d}t^2}+\left(Q_2\frac{1}{2}\frac{\mathrm{d}Q_1}{\mathrm{d}t}\frac{1}{4}Q_1^2\right)w=Q_0(t)\mathrm{exp}\left(\frac{1}{2}_0^tQ_1(\tau )d\tau \right).$$
(A3)
The study of this equation is easier than (A1).

We assume $`P^{}`$ as the dependent variable. This puts eq. (11) in the form of (A1). We then apply eq. (A2) and obtain eqs (12) and (13).


text:
# Bi2Sr2CaCu2O8+δ Bicrystal 𝑐-Axis Twist Josephson Junctions: A New Phase-Sensitive Test of Order Parameter Symmetry

## I Introduction

Phase-sensitive experiments to test the symmetry of the superconducting order parameter (OP) in the high transition temperature ($`T_c`$) superconductors were mostly made on YBa2Cu3O7-δ (YBCO) , for which the OP can have mixed ($`d_{x^2y^2}+s`$) symmetry. To reconcile the various results, Müller proposed that the surface might be mostly $`d_{x^2y^2}`$-wave, and the bulk mostly $`s`$-wave . Especially in Bi2Sr2CaCu2O8+δ (BSCCO), the $`c`$-axis transport above $`T_c`$ is incoherent, , and scanning tunneling microscope studies revealed that it is electronically disordered on the scale of the superconducting coherence length, $``$1.5 nm , both unfavorable features for bulk $`d`$-wave superconductivity. Both $`c`$-axis Pb/BSCCO Josephson junctions  and a new phase-sensitive experiment on BSCCO are consistent with those observations .

## II The Bicrystal Twist Experiment

Li et al. cleaved a single crystal of BSCCO, twisted the two cleaves an angle $`\varphi _0`$ about the $`c`$-axis and fused them together, forming junctions of remarkably superior quality . High resolution transmission electron spectroscopy and other studies revealed that the junction cross-sections were atomically clean over more than 5 $`\mu `$m , far superior to those used in tricrystal experiments . They measured the $`c`$-axis critical current $`I_c^S`$ and $`I_c^J`$ across a single crystal and the twist junction near to $`T_c`$, respectively, and the respective areas $`A^S`$ and $`A^J`$. They found that the ratio of the critical current densities $`J_c^J=I_c^J/A^J`$ to $`J_c^S=I_c^S/A^S`$ at $`0.9T_c`$ was independent of $`\varphi _0`$, as shown in Fig. 1 . Here we argue that these data demonstrate that the bulk OP in BSCCO is $`s`$-wave for $`TT_c`$, and the $`c`$-axis quasiparticle tunneling is strongly incoherent.

## III Group Theory and the Fermi Surface

Both YBCO and BSCCO are orthorhombic, but in different ways. For YBCO with point group $`C_{2v}^1`$, the mirror planes $`\sigma _x,\sigma _y`$ (the $`ac,bc`$ planes) contain the crystal axes $`a`$ and $`b`$ along the Cu-O bond direction in the CuO2 layers. In BSCCO with approximate point group $`C_{2v}^{13}`$, the mirror plane $`\sigma _b`$ (the $`bc`$-plane) contains the $`b`$ crystal axis (along a diagonal between the Cu-O bond directions) and the periodic lattice distortion . The irreducible representations for the OPs in YBCO and BSCCO are given in Table I. Although $`s`$\- and $`d_{x^2y^2}`$-wave OP components are compatible in the bulk of YBCO, they are incompatible in BSCCO, requiring a second phase transition for bulk coexistence.

We assume the quasiparticle dispersion has either the tight-binding $`\xi (𝐤)=t[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)]+t^{}\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)\mu `$ with $`t=306`$ meV, $`t^{}/t=0.90`$, and $`\mu /t=0.675`$, or hot spot $`[\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)]^2\nu ^2`$ forms, and a respective tetragonal Fermi surface (FS) with $`\xi (𝐤_F)=0`$, shown in Fig. 2 .

## IV The Twist Theorem and its Consequences

For weak tunneling across the twist junction,

$`J_c^J(\varphi _0)`$ $`=`$ $`|4eT{\displaystyle \underset{\omega }{}}f_{𝐤,𝐤^{}}^JF_{\omega ,𝐤}[R_𝐤^{}(\varphi _0)F_{\omega ,𝐤^{}}]|,`$ (1)
where $`F_{\omega ,𝐤}=\mathrm{\Delta }(𝐤)/[\omega ^2+\xi ^2(𝐤)+|\mathrm{\Delta }(𝐤)|^2]`$, $`\mathrm{\Delta }(𝐤)`$ is the OP, $`\omega `$ represents the Matsubara frequencies, $`f_{𝐤,𝐤^{}}^J`$ is the tunneling matrix element squared, $`\mathrm{}`$ is an average over each first Brillouin zone (BZ), and $`R_𝐤^{}(\varphi _0)`$ rotates the wave vectors $`𝐤^{}`$ by $`\varphi _0`$ about the $`c`$-axis. For a single $`d`$-wave OP component, $`R_𝐤(\pi /2)\mathrm{\Delta }(𝐤)=\mathrm{\Delta }(𝐤)`$.

###### Twist Theorem 1

For any weak tunneling matrix element squared satisfying $`f_{𝐤,𝐤^{}}^J=f_{𝐤^{},𝐤}^J`$, an arbitrary OP of general $`d_{x^2y^2}`$\- or $`d_{xy}`$-wave symmetry in a tetragonal crystal gives rise to a vanishing $`c`$-axis critical current across an internal 45 twist junction for $`TT_c`$.

PROOF

$`Z_\omega `$ $`=`$ $`f_{𝐤,𝐤^{}}^J[R_𝐤^{}(\pi /4)F_{\omega ,𝐤^{}}]F_{\omega ,𝐤}`$ (2)
$`=`$ $`f_{𝐤,𝐤^{}}^JF_{\omega ,𝐤^{}}[R_𝐤(\pi /4)F_{\omega ,𝐤}]`$ (3)
$`=`$ $`f_{𝐤,𝐤^{}}^JF_{\omega ,𝐤^{}}[R_𝐤(\pi /4)F_{\omega ,𝐤}]`$ (4)
$`=`$ $`f_{𝐤^{},𝐤}^JF_{\omega ,𝐤}[R_𝐤^{}(\pi /4)F_{\omega ,𝐤^{}}]=Z_\omega =0.`$ (5)
We studied a variety of OP and weak $`f_{𝐤,𝐤^{}}^J`$ forms . In Fig. 3, we show our results for the ordinary-$`d`$\- and ordinary-$`s`$\- OP forms proportional to $`\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)`$ and $`1`$, respectively. For a Gaussian form of $`f^J`$, $`f_{𝐤,𝐤^{}}^J=f_0^J\mathrm{exp}[(𝐤𝐤^{})^2/\stackrel{~}{\sigma }^2]`$. We also studied exponential, Lorentzian, rotationally-invariant Lorentzian, and stretched Lorentzian forms . Regardless of the form of $`f^J`$, the twist theorem requires $`J_c^J(\pi /4)=0`$ for a $`d`$-wave OP. From these and similar unpictured curves, it is evident that only an OP of general $`s`$-wave symmetry can fit the data. Moreover, the $`c`$-axis quasiparticle tunneling must be very incoherent.

### A Orthorhombicity

When orthorhombicity is included, the theorem is not rigorous. However, including a small $`g_{xy(x^2y^2)}`$ OP component in a dominant-$`d_{x^2y^2}`$ $`A_2^{}`$ OP will only shift the angle $`\varphi _0^{}`$ at which $`J_c^J(\varphi _0^{})=0`$ by a small amount from 45, as pictured in Fig. 4. Hence, orthorhombicity cannot explain the data of Li et al. .

### B Order Parameter Twisting

At low $`T`$, it is possible to obtain $`J_c^J(\pi /4)0`$ with a predominant-$`d_{x^2y^2}`$ $`A_2^{}`$ OP symmetry, provided that a subdominant $`A_1^{}`$ OP component can exist. Near to the twist junction, the dominant $`A_2^{}`$ OP would be suppressed, and the subdominant $`A_1^{}`$ OP increases in amplitude, so that the overall OP effectively rotates near the twist junction . However, the amount of twisting is strongly limited by the second bare transition temperature $`T_{cB}^0`$ and by the bulk and twist junction Josephson coupling strengths $`\eta ,\eta ^{}`$. In Fig. 5, we show the results obtained with subdominant OPs of the $`d_{xy}`$\- and $`s`$-wave forms, respectively. Neither case can fit the data of Li et al. . Figure 5 shows that an experiment just below $`T_c`$ can rule out OP twisting effects, unless $`T_{cB}^0=T_c`$, for which the overall bulk OP would be nearly isotropic at low $`T`$.

### C Strong Coupling and Other Models

We also considered both weak and strong coherent $`c`$-axis tunneling with a variety of FS forms . In Fig. 6, we show the tight-binding FS $`J_c^J(\varphi _0)`$ for second order coherent tunneling. In Fig. 7, the $`J_c^J(\varphi _0)`$ for weak coherent tunneling with two hot spot FSs in Fig. 2 are shown. None of these curves fit the data of Li et al..

## V BSCCO Whisker $`c`$-Axis Twist Experiments

Recently, Takano et al. performed low-$`T`$ $`c`$-axis twist mesa experiments using overdoped BSCCO whiskers with $`45^{}\alpha \varphi _090^{}`$ . Their data, pictured in Fig. 7, are distinctly different from those of Li et al. , with a strong $`J_c^J(\varphi _0)`$ dependence. Especially since $`J_c^J(\varphi _0)`$ for $`\varphi >80^{}`$ was anomalously large, Takano et al. suspected an extrinsic $`\varphi _0`$ dependence to $`f_{𝐤,𝐤^{}}^J`$ . Nevertheless, in Fig. 7 we fit the data using Eq. (1) by assuming the quasiparticles have a hot spot dispersion and intrinsically coherent $`c`$-axis tunneling. Subsequently, they found $`J_c^J(\varphi _0)C0`$ from many junctions with $`\varphi _045^{}`$, and provided preliminary Fraunhofer and Shapiro evidence that the non-vanishing $`J_c^J(45^{})`$ arises from first-order Josephson tunneling . Hence, the whisker experiments rule out a pure $`A_2^{}`$ (e. g., $`d_{x^2y^2}`$) OP, but are presently consistent with an OP either of pure $`A_1^{}`$ (e. g., $`s`$) symmetry, or of mixed $`A_1^{}`$ and $`A_2^{}`$ (e. g., $`d_{x^2y^2}+is`$) symmetry. From our theoretical studies, whisker experiments just below $`T_c`$ might determine if their bulk OP is also pure $`A_1^{}`$ , and measurements above $`T_c`$ could investigate if the $`c`$-axis transport is indeed coherent, strikingly inconsistent with single crystal BSCCO .

## VI Conclusions

The data of Li et al. demonstrate that the OP in the bulk of BSCCO has $`A_1^{}`$ ($`s`$) symmetry for $`TT_c`$, and that the $`c`$-axis tunneling is strongly incoherent. The data of Takano et al. presently rule out a pure $`A_2^{}`$ ($`d_{x^2y^2}`$) OP, but surprisingly suggest that the weak $`c`$-axis tunneling in BSCCO whiskers might be coherent with a FS consisting of small hot spots.

## acknowledgments

We thank G. B. Arnold, Qiang Li, K. A. Müller, K. Scharnberg, and M. Tachiki for useful discussions.