text:
# A Kinematic Model for the Narrow-Line Region in NGC 4151Based on observations made with the NASA/ESA Hubble Space Telescope. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS5-26555.

## 1 Introduction

The study of the kinematics of the narrow line region (NLR) in Seyfert galaxies is important, because it has the potential to reveal the dynamical forces at work on the emission-line gas within about one kpc of the nucleus. In this region, the central supermassive black hole presumably dominates the kinematics, due to its gravitational influence and/or the radiation, winds, and jets that emanate from its immediate surroundings. With the Hubble Space Telescope (HST), we have the ability to resolve the emission-line knots within the NLR and with the two-dimensional detectors of the Space Telescope Imaging Spectrograph (STIS), we can measure their velocities accurately and efficiently. In this paper, we concentrate on STIS spectra of the apparently brightest Seyfert 1 galaxy, NGC 4151, and compare our results to those from a kinematic study based on STIS spectra of the brightest Seyfert 2 galaxy, NGC 1068 (Crenshaw & Kraemer 2000a, hereafter CK2000).

The NLR of NGC 4151 as defined by HST \[O III\] images (Evans et al. 1993) is biconical in shape with the apex located at the bright nucleus, and contains numerous bright emission-line knots and filaments. Evans et al. find a projected half-opening angle of 37.5o $`\pm `$ 5o for the bicone and a position angle (PA) of the bicone axis on the sky of 60o $`\pm `$ 5o. The extended narrow-line region (ENLR), at distances greater than 6′′ from the nucleus, lies within the extrapolated bicone, but most of the emission from this region is concentrated between position angles 30o and 60o (Pérez-Fournon & Wilson 1990). Ground-based spectroscopy has revealed that the \[O III\] and H I velocity fields are very similar at these distances (Pedlar et al. 1992), indicating that the kinematics in the ENLR are dominated by galactic rotation.

Although the velocity field of the NLR in NGC 4151 (at distances less than 6′′ from the nucleus) is difficult to determine from the ground, Schulz (1990) found evidence for outflow along the cones, superimposed on galactic rotation. With observations primarily from HST’s Faint Object Camera (FOC), which has a limited ability to obtain long-slit spectra, Winge et al. (1997, 1999) claimed that the dominant kinematic component is rotation. STIS observations and detailed kinematic models are crucial for testing these claims. As part of a STIS Key Project carried out by members of the Instrument Definition Team, we have accumulated a number of observations of NGC 4151. Our previous results on these data are described in Hutchings et al. (1998, 1999), Kaiser et al. (2000), Nelson et al. (2000), and Kraemer et al. (2000). We describe the observations that are pertinent to our kinematic model in the next section.

We take the systemic velocity of NGC 4151 to be 997 km s-1, based on H I observations of the outer regions of the host galaxy (Pedlar et al. 1992). For a Hubble constant of H0 $`=`$ 75 km s-1 Mpc-1 and a distance to NGC 4151 of 13.3 Mpc, 0$`^{\prime \prime }.`$1 corresponds to 6.4 pc.

## 2 Previous STIS Observations and Results

Our first STIS observation was a slitless spectrum of the H$`\beta `$ and \[O III\]$`\lambda \lambda `$4959, 5007 emission at a spectral resolving power of $`\lambda `$/$`\mathrm{\Delta }\lambda `$ $``$ 8000 (Hutchings et al. 1998). By comparing the slitless spectrum with an archival WFPC2 image obtained through a narrow-band \[O III\] filter, we were able to determine the radial velocities and velocity dispersions of about 40 bright emission-line knots (with sizes on the order of 0$`^{\prime \prime }.`$2). We could also identify about 10 knots in a slitless spectrum of the C IV region, and found that their velocities were consistent with the \[O III\] measurements. We found that the highest velocities are close to the nucleus, and, in general, the emission-line knots are blueshifted to the southwest and redshifted to the northeast. A rotation curve cannot match the observed radial velocities as a function of distance. We proposed a simple model consistent with the geometry of Evans et al. (1993) to explain the gross kinematics: radial outflow along the surface of a bicone, with an opening angle of 80o and a line of sight that is 10o outside of the cones.

In Kaiser et al. (2000), we presented slitless spectra of several spectral regions at a more favorable roll angle than that of the Hutchings et al. (1998) observations. Using these data, we were able to determine radial velocities for a number of additional emission-line knots. In Figure 1, we compare the radial velocities for the emission-line knots in common in these two studies. The average difference in radial velocity is only 47 km s-1, which indicates the reliability of this observational technique. In Kaiser et al. (2000), an attempt to fit the radial velocities as a function of distance with a rotating disk model and point-source mass was not satisfactory; the transition between blueshifted and redshifted velocities is much steeper than observed, and an unreasonably large central mass is required: 1010M if all clouds are included, and 4 x 108M if clouds with radial velocities of only 200 km s-1 or less are included. Thus, we found that outflow is extremely likely, since gravitational motions (including infall) cannot account for the high velocities. There is no strong correlation between velocity or velocity dispersion and the positions of the radio knots, which indicates that jet acceleration is not important; wind or radiation pressure are more likely sources of the acceleration.

In Hutchings et al. (1999), we used narrow-band WFPC2 images of NGC 4151 to isolate faint line emission at high velocities for several ions. Radial velocities up to 1400 km s-1 were seen in both approach and recession on both sides of the nucleus, indicating a separate kinematic component from that seen in the bright emission-line knots. We will not attempt to model this component in the current paper; additional progress on this component is likely to come from planned long-slit observations at high spectral resolution.

In Nelson et al. (2000), we presented long-slit (52′′ x 0$`^{\prime \prime }.`$1) spectra of NGC 4151, obtained at low spectral resolving power ($`\lambda `$/$`\mathrm{\Delta }\lambda `$ $``$ 1000) over the 1150 – 10,000 Å range at two position angles: 221o (passing through the nucleus) and 70o (passing 0$`^{\prime \prime }.`$1 south of the nucleus, see Figure 1 in Nelson et al. 2000 for the slit placements superposed on WFPC2 \[O III\] images). Although these data do not have the spatial coverage and spectral resolution of the slitless spectra, the slit allowed us to determine the location and radial velocities of much weaker emission, including that between the bright emission-line knots. To determine the radial velocities, we used the \[O III\] $`\lambda `$5007 emission in 0$`^{\prime \prime }.`$2 bins along the 0$`^{\prime \prime }.`$1-wide slit. We explored two basic models in this paper: radial outflow away from the nucleus and expansion away from the radio jet axis. We concluded that radial outflow provides the best match with the data.

## 3 Data Analysis and Measurements

We rely on our previous observations to develop a detailed kinematic model for the NLR of NGC 4151 in this paper. For our analysis of the slitless spectra, we use the measurements of Kaiser et al. (2000), which we have shown are consistent with those of Hutchings et al. (1998). For the long-slit spectra, we remeasured the \[O III\] $`\lambda `$5007 emission along the slit at each pixel location (0$`^{\prime \prime }.`$05 intervals), rather than in 0$`^{\prime \prime }.`$2 bins, to avoid smoothing of features in the radial velocity curves. At some locations close to the nucleus, there are two distinct components of the \[O III\] emission at this spectral resolution. We fit the spectrum at each location with a local continuum plus a Gaussian for each clearly identifiable peak of \[O III\] emission, resulting in one or two kinematic components. We determined uncertainties in the radial velocities of the components from the errors in the Gaussian fits and from different reasonable continuum placements.

Figures 2 and 3 give the radial velocities, widths (full-width at half-maximum, corrected for the instrumental value of 300 km s-1), and fluxes as a function of angular distance from the nucleus for the long-slit spectra at the two position angles; our convention is that negative positions represent the southwest direction. The radial velocity curves in the top plots are very similar to those determined by Nelson et al. (2000), except they show a little more structure due to the finer sampling. At distances less than 4′′, the emission-line knots southwest of the nucleus are blueshifted and the knots northeast of the nucleus are redshifted. At greater distances, the radial velocities approach the rotational velocities seen in the ENLR. High radial velocities (greater than $``$300 km s-1) tend to occur within $``$1′′ of the nucleus, although low radial velocities occur in this region as well. In many cases, the structure in the radial velocity curves can be attributed to the varying contributions of emission-line knots and background emission to the total \[O III\] emission as a function of position. The knotty structure can be seen in the bottom plots of Figures 2 and 3. The middle plots in Figures 2 and 3 show that there are no strong trends in the velocity widths, except that the lower velocity clouds at large distances tend to have smaller widths. The velocity widths are due to superposition of emission knots with distinct velocities as well as internal dispersion within each knot.

In Figure 4, we compare the radial velocities derived from the long-slit spectra with those from the slitless spectra; the number of points are limited, since there only a few bright emission-line knots that lie along the slit locations. Again, the measured radial velocities are very similar, with no evidence for systematic offsets as a function of velocity or dispersion. The average difference in radial velocity is only 37 km s-1. We conclude that the two observational techniques provide reliable radial velocities; slitless spectra at moderate resolution are useful for obtaining velocities of all of the bright emission-line knots in a single observation, whereas long-slit low-resolution spectra are helpful for determining the kinematics of all of the emission along a particular direction.

## 4 A Biconical Outflow Model

Our model for NGC 4151 is based on a previous kinematic model for the NLR in the Seyfert 2 galaxy NGC 1068 (CK2000). From STIS low-resolution long-slit spectra of the \[O III\] emission, we found that a simple model of biconial radial outflow provides a good match to the radial velocities as a function of distance from the hot spot in NGC 1068. In that model, the bicone is evacuated along its axis, which is inclined only 5o out of the plane of the sky. The emission-line knots show radial acceleration to a distance of $``$140 pc from the hot spot, followed by deceleration that approaches the systemic velocity at $``$310 pc. The deceleration is likely due to collision with a patchy and anisotropically distributed ambient medium.

For the biconical model of NGC 4151, we assume that the two cones have identical properties (geometry, velocity law, etc.), a filling factor of one within the defined geometry, and no absorption of \[O III\] photons. The parameters that are allowed to vary in our code are the extent of each cone along its axis (zmax), its minimum and maximum half-opening angles ($`\theta `$inner and $`\theta `$outer), the inclination of its axis out of the plane of the sky (iaxis), and the velocity law as a function of distance from the nucleus. Our code generates a two-dimensional velocity map and samples this map with a slit that matches the position, orientation, and width of the observational slits. Note that we are only attempting to match the observed radial velocities with the model, and not the observed line widths or fluxes.

Given the success of our model for NGC 1068, we assumed the same velocity law for NGC 4151: constant acceleration to a maximum velocity (vmax) at a turnover radius (rt), followed by a constant deceleration to zero velocity at the greatest extent of the cone ($`=`$ zmax/cos $`\theta `$outer). Then we adjusted the other parameters of the model to determine if a reasonable match to the observed radial velocities could be obtained. First, we tilted the bicone axis out of the sky until the observed form of the model radial velocity curves matched that of the observations. For NGC 4151, this occurs when the bicone moves completely out of the plane of the sky, so that one cone is blueshifted and the other side is redshifted. Then, we varied the other parameters until we obtained a reasonable match with the observations. Since the observed radial velocity curves show significant scatter, we did not attempt to fine-tune the models; instead we settled for an illustrative model that matches the overall trends.

The parameters for our kinematic model of NGC 4151 are shown in Table 1; values for the model of NGC 1068 are given for comparison. The bicone opening angles, maximum extent, and turnover radius are very similar for the two models. The maximum velocity at the turnover is much lower in NGC 4151, which may be due to its lower intrinsic luminosity (see the Discussion). The principal difference between the two models is the viewing angle, with the bicone axis of NGC 4151 much closer to the line of sight. However, as we note in the Discussion, the inclination is not enough to place the line of sight into the ionization cone. Given the inclination angle, we find that the projected maximum half-opening angle onto the plane of the sky is 43.5o, which is close to the 37.5o$`\pm `$5o determined from the HST images by Evans et al. (1993).

Two-dimensional velocity maps for the model of NGC 4151 are given in Figure 5, along with the slit positions and widths that correspond to our long-slit observations. For demonstration, we show only the front and back sides of the outermost surface (at $`\theta `$outer). For comparison, we show the same surfaces in Figure 6 for our model of NGC 1068 (not shown in CK2000 due to space limitations). In both figures, the models are aligned so that the projected bicone axis on the sky is vertical, and northeast is at the top. The inclination of the bicone axis is NGC 4151 is such that the entire northeast cone is redshifted, and the entire southwest cone is blueshifted. By contrast, the axis for NGC 1068 is inclined by only 5o so that blueshifts and redshifts of similar magnitude occur in both cones (see CK2000 for more details). The effects of our acceleration plus deceleration velocity laws, modified by projection of the velocity vectors from the bicone, are visible in these images. One interesting feature of the NGC 4151 model is that the maximum radial velocities occur very close to the nucleus on the plane of the sky, regardless of slit position.

## 5 Comparison of the Observations and Models

First, we compare the models with the long-slit observations. Figure 7 shows the envelopes of radial velocities from the model of NGC 4151, extracted from the two slit positions, compared to the observed radial velocities (relative to systemic). There are four distinct components of the envelope in each plot, corresponding to the four quadrants of the bicone intersected by the slit. The width of the envelope is determined by the range in half-opening angle, and the relative amplitudes are determined by the inclination of the bicone axis. Although the match is not perfect, the overall trends are reproduced. The major discrepancy is that this model does not explain the blueshifted velocities in the $``$150 to $``$300 km s-1 range (relative to systemic) in the southwest cone. This discrepancy could be eliminated by filling in the southwest cone with emitting material, which would fill in the gap in the plots, but would also fill in regions that are not occupied by observed points. Thus, the bicone for the most part must be evacuated along its axis . Our model also does not explain the emission knots with very high velocities ($`>`$ 700 km s-1 relative to systemic), which is also the case for NGC 1068. As noted in CK2000, this indicates that some emission-line knots do not experience a significant deceleration.

Next, we compare the models with the slitless observations. In the top plot of Figure 8, we show the slitless radial velocities as a function of position, relative to the nucleus at the origin. The magnitude of the radial velocity is represented by the size of the symbol, with “+” representing blueshift and “x” representing redshift, relative to the systemic velocity. The plot is oriented such that the bicone axis, at PA = 60o, is vertical. The trends discussed earlier can be seen in this plot (for a color-coded velocity map, see Hutchings et al. 1998).

The bottom plot in Figure 8 shows the observed radial velocities minus the model values. At each observed point, we have two model values to choose from: the average radial velocity for the front side and the average value for the back side; we adopt the value closest to the observed radial velocity and subtract. Although our kinematic model was calculated to fit the long-slit data, the small residual velocities in this plot show that it provides a reasonable fit to the slitless data as well. There is a mixture of positive and negative residuals in each cone, indicating the overall agreement between the model and observations. The two largest residuals in this plot identify the “rogue” blueshifted knots that apparently do not experience any deceleration, which correspond to knots 23 and 26 in Hutchings et al. (1998) and Kaiser et al. (2000).

In Figure 9 we plot the observed radial velocities vs. the model values at each point to show the magnitude of the residuals. Except for a few points, notably the two highly blueshifted knots, the residuals are relatively small and evenly distributed around the unity line. We conclude that our kinematic model provides a good match to the general trends seen in all of our slitless and long-slit data sets.

## 6 Discussion

The details of our kinematic model of NGC 4151 were inspired by STIS observations and a kinematic model of the NLR in NGC 1068. Since the outflow axis is nearly in the plane of the sky in NGC 1068, the trend of increasing velocity with distance to a turnover radius, followed by decreasing velocity, can be easily seen in this object, and the lack of emission along the cone axis in NGC 1068 is also clear (CK2000). The detailed kinematics of the NLR in NGC 4151 have been more difficult to interpret, since it turns out that the cone axis is significantly out of the plane of the sky, and velocities from different radial distances are projected along the same line of sight (see Figure 5). However, by adopting the basic geometry and velocity law used for NGC 1068, we are able to construct a kinematic model of the NLR in NGC 4151 that matches the observed trends in radial velocity as a function of position. The only major changes were to reduce the maximum radial velocity from 1300 to 750 km s-1, and to tilt the bicone axis from 5o to 40o out of the plane of the sky. Thus, we have a basic model that can explain the NLR kinematics in both Seyfert galaxies.

### 6.1 Comparison with Other Models

Winge et al. (1997, 1999) claim that the dominant kinematic component in the NLR of NGC 4151 is rotation. However, there are several problems with this model. First, as we discussed earlier, the observed transition between blueshifted and redshifted radial velocities is much shallower than the predicted values from a Keplerian disk with a point-source mass. Thus, in addition to a point-source mass of 5 x 107 M$``$, these authors must include an extended mass of order 109 M in the inner 0$`^{\prime \prime }.`$1, to match the shallow transition between blueshifted and redshifted velocities. To date, there is no confirmatory evidence for such an extended mass. Second, Winge et al. (1999) are only able to fit their low-velocity component (within 300 km s-1 of systemic) with this model; the higher velocities require another explanation. Finally, we note that a rotation model cannot explain the kinematics of the NLR in NGC 1068, which shows blueshifts and redshifts of similar magnitude on both sides of the nucleus.

To explore the uniqueness of our model for NGC 4151, we have experimented with a number of different velocity laws and inclination angles, but with the same biconical geometry. We summarize our findings as follows. 1) Simple radial outflow at constant velocity or acceleration cannot explain the decreasing velocities at distances $`>`$ 1′′. 2) Gravitational infall cannot produce the small radial velocities seen within 1′′ of the nucleus. Also, as we pointed out in Kaiser et al. (2000), gravitational infall models require unreasonably high masses (on the order of 1010 M$``$) to produce the high velocities that are seen. 3) Axon et al. (1998) propose a model for NGC 1068 in which the gas expands away from the radio jet (which is nearly coincident with the ionization cone axis). The simplest version of this model is one in which the velocity vectors are perpendicular to the bicone axis. As demonstrated by Nelson et al. (2000), this model cannot work for NGC 4151. As the bicone axis is tilted out of the plane of the sky, the velocity vectors in each of the two cones decrease in magnitude by the same amount, and there is no inclination angle that will yield the observed blueshifts on one side of the nucleus and redshifts on the other. 4) Nelson et al. (2000) show that a reasonably good match to the STIS observations can be obtained by assuming radial outflow, with velocity decreasing with distance as r-1/2. We verify that this velocity law provides an acceptable solution for NGC 4151, except that there is some indication of acceleration of the clouds away from the nucleus at small distances (see Figure 7). We note that this velocity law does not work for NGC 1068, which shows clear acceleration of the clouds to an angular distance of $``$1.7′′ (CK2000). Thus, our proposed model provides the best explanation at present for both Seyfert galaxies.

### 6.2 Geometry of the Bicone and Host Galaxy

The host galaxy of NGC 4151 is nearly face-on (inclination from our line of sight $`=`$ 20o, PA of the major axis $`=`$ 22o, see Pedlar et al. 1992). Given the inclination and PA of our kinematic model of the NLR, the bicone axis is inclined by $``$40o with respect to the host galaxy, and one side of each cone (corresponding to a low radial velocity quadrant in Figure 7) lies close to the galactic plane (similar to the case for NGC 1068, see CK2000). Thus, our proposed geometry allows the gas in the galactic disk to be ionized by the continuum source and is compatible with the idea that this gas is the ENLR. The general agreement between ground-based \[O III\] and H I radial velocities at a distance of 10 – 20′′ from the nucleus provides strong evidence for this idea (Pedlar et al. 1992); in particular, these observations indicate normal galactic rotation, with typical radial velocities of $``$70 km s-1 in the southwest and $`+`$70 km s-1 in the northeast. Our slitless observations are consistent with these values; they indicate an average radial velocity of $``$50 $`\pm `$ 15 km s-1 relative to systemic for the \[O III\] emission 10 – 20′′ southwest of the nucleus (after correcting for an offset of $``$40 km s-1 for the \[O III\] emission at the nucleus, relative to the systemic, see Kaiser et al. 2000). Thus, our observations confirm that the NLR and ENLR are two physically distinct regions in NGC 4151; from the kinematics, it appears that the transition from one region to the other occurs in the region 6′′ to 10′′ from the nucleus.

The standard unified model for Seyfert galaxies posits a thick molecular torus that produces ionization cones with sharp edges (Antonucci 1993). In this scheme, our line of sight to a Seyfert 1 galaxy should be inside the ionization bicone. However,the sum of the inclination and maximum half-opening angles from our model of NGC 4151 is 76o which places our line of sight 14o outside of the bicone. This discrepancy was pointed out Evans et al. (1993), who arrive at a viewing angle that is $``$30o outside of the bicone, based on the observed morphologies of the NLR and ENLR. Thus, NGC 4151 does not appear to fit the standard unified model. However, if the obscuring torus is not optically thick at all polar angles (cf., Evans et al. 1993; Konigl and Kartje 1994), then the ionization cones would have “fuzzy” edges, and NGC 4151 would offer no problems for this type of model.

### 6.3 Implications for Dynamical Models

Based on our observations and kinematic models, we find strong evidence for radial outflow in the NLRs of NGC 4151 and NGC 1068. The velocity flow is well-organized, suggesting that the emission-line clouds originate close to the nucleus and are forced outward from there. A linear increase in velocity with distance to a turnover radius provides a reasonable fit to the observations. There are two plausible interpretations of this trend. One is that the clouds in these two NLRs originated from an explosive outburst that occurred $``$105 years ago (the crossing time to a distance of $``$150 pc), which would naturally lead to a linear increase of velocity with distance. If this model is generally correct, then the NLRs in some Seyfert galaxies should be evacuated near the nucleus; so far there is no evidence for this in HST \[O III\] images (cf., Schmitt & Kinney 1996).

Another interpretation assumes a steady-state model, in which clouds are continuously ejected from the nucleus and experience approximately constant acceleration out to the turnover radius. Under this interpretation, the most likely sources of acceleration would be radiation and/or wind pressure. Radiative driving of resonance lines is an attractive means of accelerating the clouds, since a strong radiation field is obviously present. However, there is also evidence for outflowing winds in the form of highly ionized gas seen in absorption in the UV (Crenshaw et al. 1999) and X-rays (Reynolds 1997; George et al. 1998). Since the turnover radius is similar for the two Seyfert galaxies, the lower maximum velocity in NGC 4151 indicates a smaller acceleration. If the acceleration mechanism is radiation pressure, the lower maximum velocity in NGC 4151 may be due to a much lower luminosity, which is $``$10% of the intrinsic luminosity of NGC 1068 (Pier et al. 1994). In any event, future dynamical models will have to explain the roughly constant acceleration of the NLR clouds to the turnover point.

It is unlikely that the radio jets are responsible for the radial acceleration of the emission line knots, since there is no apparent correlation of their radial velocities with radio knot location in either NGC 4151 (Kaiser et al. 2000) or NGC 1068 (CK2000). However, the radio jets may be responsible for evacuating the bicones along their axes, since the emission-line knots appear to be avoiding the radio knots in the inner NLRs of these two Seyferts (Kaiser et al. 2000; Capetti et al. 1997).

An unexpected result of this study and that of CK2000 is that the NLR clouds in NGC 4151 and NGC 1068 undergo deceleration at a distance of $``$150 pc from the nucleus. The deceleration cannot be due to gravity because an unreasonably high mass ($``$1010 M$``$) is required at the distances and velocities of the turnover point. Following CK2000, the simplest explanation for the deceleration is that the clouds experience a drag force due to a surrounding ambient medium extending outward from the turnover point. In this interpretation, the emission knots with unusually high velocities have not experienced drag forces because they have passed through holes in the surrounding medium. Additional evidence for collision of the NLR clouds with an ambient medium comes from a study of the physical conditions in the NLR of NGC 1068 (Kraemer & Crenshaw 2000). Near the location of the velocity turnover in the blueshifted region of the northeast cone, we find evidence for an additional source of EUV ionizing radiation that is likely to be generated from the shock front as the clouds encounter the ambient medium. However, we note that these effects have not (as yet) been detected at other turnover points in NGC 1068 or in NGC 4151. Another clue to the ambient medium can be found from observations of the extended continuum in NGC 1068 (Crenshaw & Kraemer 2000b), which show very enhanced regions of electron-scattered continua. As pointed out by Kraemer & Crenshaw (2000), these are likely to be regions of highly-ionized, low-density (ionization parameter $``$ 1, electron density $``$ 100 cm-3) gas which must nearly fill the bicone at these distances, and are therefore likely candidates for the decelerating media. Unfortunately, the glare from the nucleus of NGC 4151 prevents us from searching for a scattered continuum component in its NLR.

## 7 Conclusions

The principal advantage of a kinematic model that invokes radial outflow in the NLR is that it can explain the observed radial velocities in both a Seyfert 1 galaxy (NGC 4151) and a Seyfert 2 galaxy (NGC 1068). Assuming that radial outflow is common in the NLR of all Seyfert galaxies, it is interesting (but not surprising) that the kinematic signature of circumnuclear gas within $``$500 pc of an active nucleus is probably different from that surrounding “inactive” supermassive black holes in nearby galaxies, in which gravitational motions dominate (e.g., M84, Bower et al. 1998). This may indicate that once a supermassive black hole “turns on” (e.g., due to an onset of fueling), radiation and/or wind pressure from the active nucleus will dominate the kinematics of the NLR.

It is intriguing that very similar geometries and velocity laws can explain the observed radial velocities in the two Seyfert galaxies that we have studied. This may indicate a common dynamical situation in all Seyferts, or it may just be a coincidence. New STIS observations of more Seyfert galaxies will be helpful in testing the general applicability of these kinematic models and providing the refinements needed for detailed dynamical models.

This work was supported by NASA Guaranteed Time Observer funding to the STIS Science Team under NASA grant NAG 5-4103.


format: text

text:
# Statistics of stable marriages

## I Introduction

The study of complex systems that are built up of many individuals interacting according to relatively simple laws has recently attracted a lot of interest among physicists. Concepts and methods developed in statistical physics have successfully been applied to describe collective and nonlinear phenomena arising in such systems. Studies on traffic flow, stock markets, voting, flocking birds, evolution etc. are only a few examples to mention here.

In the stable marriage problem the elements of two sets have to be matched by pairs. As possible applications one could think of job seekers and employers, lodgers and landlords, or simply men and women who want to get married. For the sake of clearness and simplicity we will exclusively refer to the paradigm of marriage in this paper. The purpose of most studies of the stable marriage problem has so far been to unravel the underlying mathematical structure and to develop new or improve existing algorithms for finding stable marriages . Only recently it has been recognized that the stable marriage problem has many features in common with classical disordered systems in statistical physics like spin glasses, in particular the element of frustration and the existence of many competing states . Applying methods developed in statistical mechanics, Nieuwenhuizen  solved several variants of the stable marriage problem, including bachelors and polygamy. In his approach, Nieuwenhuizen focussed on the properties of globally optimal matchings which are advantageous for the society as a whole, but not necessarily for all individuals. Here we consider the opposite limit where all individuals pursue only their own egoistic objectives regardless of the consequences for others. Under these conditions, only matchings that are stable with respect to the actions of all individuals are relevant. In the framework of game theory this kind of stability is termed Nash equilibrium.

Our paper is organized as follows. In Section 2 we introduce the stable marriage model and set up the notation. The Gale-Shapley algorithm and the properties of the stable matchings obtained with it are presented in Section 3. In Section 4 we analyze the set of all stable matchings. Generalizations of the model, introducing an acceptance threshold and unequal number of men and women are discussed in Sections 5 and 6, respectively. Finally, we summarize our results in Section 7.

## II The stable marriage model

We consider a society consisting of $`N`$ men and $`N`$ women who have to be matched by pairs. The cost for man $`i`$ to be married with woman $`j`$ is $`x(i,j)`$ and the cost for woman $`j`$ to be married with man $`i`$ is $`y(j,i)`$. In order to investigate statistical properties of the model we assume that all cost functions are independent random variables uniformly distributed between $`0`$ and $`1`$ In Ref.  a different probability distribution has been considered, namely $`\rho (x)x^r\mathrm{exp}(x)`$. The case $`r=0`$ corresponds essentially to our choice.. Since all individuals simultaneously attempt to optimize their own benefit by being married with a partner of lowest possible cost, conflicting wishes arise, and in general it is impossible to find a matching where everybody is perfectly satisfied. Let us assume that a matching is formed where man $`i`$ is married with woman $`j`$. We define $`x_ix(i,j)`$ and $`y_jy(j,i)`$, i.e. $`x_i`$ is the cost for man $`i`$ in this particular matching and $`y_j`$ the corresponding cost for woman $`j`$. Instead of cost we will frequently use the expression energy in order to emphasize the analogy with physical systems where the energy is the quantity to be minimized in the groundstate. For later convenience, we introduce the total energies of men and women by

$$X=\underset{i=1}{\overset{N}{}}x_iY=\underset{j=1}{\overset{N}{}}y_j$$
(1)
In search of the rules according to which marriages can be formed or broken, essentially two points of view have emerged: In the first one, the goal is to find the globally best matching, i.e. the one where the total energy $`U=X+Y`$ is minimal. This problem has recently been studied by Nieuwenhuizen  using techniques borrowed from statistical physics of disordered systems, such as the replica trick. The second point of view, that we will adopt, emphasizes the role of individuals as decision-makers. As long as there is no supervising body like church or government that dictates who has to be married with whom no individual is obliged to stay in an unsatisfactory relationship if he or she can find a better one. This leads us to the notion of stability. A matching is called stable if there is no man and no woman who prefer each other to their actual partners. In other words, a matching is unstable, if there exists a man $`i`$ and a woman $`j`$, who are not married with each other, such that

$$x_i>x(i,j)\mathrm{and}y_j>y(j,i)$$
(2)
The particular significance of stable matchings lies in the fact that they cannot be broken, neither by the action of a single pair nor by any coalition of men and women. In contrast to the globally best matching which can only be found with the help of computers for very small numbers of men and women, there exist powerful algorithms to search for stable matchings, even for systems of several thousands of individuals. The simplest of these algorithms, introduced by Gale and Shapley some time ago , will be presented in the next section.

## III The Gale-Shapley algorithm

In the man-oriented Gale-Shapley (GS) algorithm only men make proposals while women accept or reject. In the beginning every man makes up a list where he places the women according to his preferences. The smaller the cost $`x(i,j)`$ the higher the position of woman $`j`$ on the list of man $`i`$. Then a series of proposals starts. Whenever a man is not married he makes a proposal to the woman with the highest rank among those to whom he has not yet proposed. If the woman is not engaged she accepts the proposal and they get married. If the woman is already married, she accepts if she prefers the suitor to her husband. In this case she divorces and gets married with the proposer. The algorithm stops when the last woman has received an offer and gets married. The number of proposals $`r`$ that a man has to make on the average during the execution of the GS algorithm has been calculated in . The result is

$$r\mathrm{log}N+\mathrm{C}$$
(3)
where $`\mathrm{C}=0.5772\mathrm{}`$ is Euler’s constant and corrections are negative and of order $`(\mathrm{log}N)^2/N`$. Amazingly, the knowledge of $`r`$ is already sufficient to calculate the total energy of men ($`X`$) and women ($`Y`$) in the GS matching. With every proposal the energy of one man increases by $`1/N`$ since he has to move forward one step on his preference list. Consequently, the total energy of men is on the average the same as the number of proposals that every man makes, i.e. $`X=r`$.

Each woman also receives on the average $`r`$ proposals. Keeping only the best offer i.e. the smallest out of $`r`$ random numbers between $`0`$ and $`1`$ yields a best value of $`y=1/(r+1)`$. Taking into account that the number of proposals that a woman receives is not fixed, but distributed according to a binomial distribution yields $`y=(1\mathrm{exp}(r))/r1/r`$ (see Appendix B). Using Eq. (3) we obtain the following estimates for the energies of men and women in the man-oriented GS matching

$$X\mathrm{log}N+\mathrm{C}Y\frac{N}{\mathrm{log}N+\mathrm{C}}$$
(4)
which implies the relation $`XY=N`$.

Fig. 1 shows the average energies $`X`$ and $`Y`$ of men and women obtained with the man-oriented GS algorithm as function of $`N`$. Each data point represents an average over several thousand realizations. A logarithmic scale is used for the energy axis in order to highlight the good agreement with the analytical result of Eq.(4) even for small values of $`N`$.

Of course, everything that has been said about the man-oriented GS algorithm is also valid in the woman-oriented version where the roles of men and women are interchanged, i.e. women propose and men judge. It has been shown  that the stable matching obtained from the man–oriented execution of the GS algorithm is man–optimal in the sense that no man can have lower energy in any other stable matching. Similarly, it is the worst possible stable matching for women. Therefore the energies obtained by performing the GS algorithm provide an upper and a lower bound for the possible energy range of all stable matchings. Furthermore, if the man and the woman–optimal matching happen to be identical we know for sure that in this case there exists only one single stable matching.

## IV The set of all stable matchings

If $`N`$ is large, there are generally many other stable matchings besides the man and the woman-optimal matchings that can be obtained with the GS algorithm. Although the maximum number of stable matchings in a system of size $`N`$ is not known, there exist lower bounds ; e.g. for $`N=32`$ the maximum number of stable matchings is larger than $`10^{11}`$ (!).

In order to calculate the average number of stable matchings we consider an arbitrary matching of $`N`$ men and $`N`$ women. This matching is unstable, if there exists a man $`i`$ and a woman $`j`$ who are not married with each other, and whose mutual costs $`x(i,j)`$ and $`y(j,i)`$ are both lower than the ones in the existing matching, $`x_i`$ and $`y_j`$, respectively. The probability that man $`i`$ and woman $`j`$ prefer to stay with their respective partners is therefore $`p_{ij}=1x_iy_j`$ and the whole matching is stable with probability

$$P=_0^1\mathrm{d}^Nx_0^1\mathrm{d}^Ny\underset{ij}{}(1x_iy_j)$$
(5)
Although this $`2N`$-fold integral, originally derived by Knuth , looks quite simple, it cannot be done exactly unless $`N`$ is very small. Using rather involved probabilistic arguments, Pittel  has derived an asymptotic formula,

$$P\frac{\mathrm{log}N}{\mathrm{e}\mathrm{\Gamma }(N)}$$
(6)
valid in the limit $`N\mathrm{}`$. However, comparing Eq. (6) with numerical simulations we found rather serious discrepancies (see Fig. 2). In order to remove these discrepancies and to provide a more transparent approach we replace the factors $`1x_iy_j`$ in Eq. (5) by $`\mathrm{exp}(x_iy_j)`$, which is justified since the dominant contribution to the integral comes from the region where the products $`x_iy_j`$ are small. Dropping the constraint $`ij`$ we arrive at the approximation

$$\stackrel{~}{P}=_0^1\mathrm{d}^Nx_0^1\mathrm{d}^Ny\mathrm{exp}(XY)$$
(7)
where the variables $`X=_{i=1}^Nx_i`$ and $`Y=_{j=1}^Ny_j`$ are the energies of men and women, respectively. Eq. (7) can also be viewed as the partition sum of two species of particles confined to the interval $`[0,1]`$ whose interaction is the product of their center of mass coordinates. It is convenient to replace the integral over the variables $`x_i`$ and $`y_j`$ by an integral over $`X`$ and $`Y`$

$$\stackrel{~}{P}=_0^NdX_0^NdY\rho (X)\rho (Y)\mathrm{exp}(XY)$$
(8)
where the probability distribution of $`X`$ is asymptotic to

$$\rho (X)\frac{X^{N1}}{\mathrm{\Gamma }(N)}\left(1\mathrm{exp}(\frac{N}{X})\right)^N$$
(9)
for small $`X`$ (see Appendix A). Inserting Eq. (9) in Eq. (8) and introducing $`t=XY`$ as new integration variable we obtain

$`\stackrel{~}{P}`$ $`=`$ $`\mathrm{\Gamma }(N)^2{\displaystyle _0^{\frac{N}{\mathrm{log}N}}}dXX^1{\displaystyle _0^{\frac{NX}{\mathrm{log}N}}}dtt^{N1}\mathrm{exp}(t)`$ (10)
$``$ $`\mathrm{\Gamma }(N)^1{\displaystyle _{\mathrm{log}N}^{\frac{N}{\mathrm{log}N}}}dXX^1`$ (11)
$`=`$ $`\mathrm{\Gamma }(N)^1\left(\mathrm{log}N2\mathrm{log}(\mathrm{log}N)\right)`$ (12)
where we have assumed that the factor $`\left(1\mathrm{exp}(\frac{N}{X})\right)^N`$ imposes a cut-off $`N/\mathrm{log}N`$ on the $`X`$-integration and the same holds for $`Y`$. The integral over $`t`$ yields $`\mathrm{\Gamma }(N)`$, provided that $`X>\mathrm{log}N`$. Comparison with Eq. (6) shows that the ratio $`\stackrel{~}{P}/P`$ approaches $`\mathrm{e}=\mathrm{exp}(1)`$ in the limit $`N\mathrm{}`$ which we also checked numerically. Since there are in total $`N!=\mathrm{\Gamma }(N+1)`$ matchings the average number of stable matchings is given by

$$S\frac{N}{\mathrm{e}}\left(\mathrm{log}N2\mathrm{log}(\mathrm{log}N)\right)$$
(13)
$`S`$ is plotted as function of $`N`$ in Fig. 2. The circles are data obtained from numerical simulations using the algorithm of Gusfield and Irving  that allows the determination of all stable matchings. Each data point represents an average over 200 random realizations of the preference lists. The dashed curve is the result of Eq. (13) and the solid curve is the asymptotic formula derived by Pittel . Note the large error made by neglecting the $`\mathrm{log}(\mathrm{log}N)`$ term.

In the following we investigate how the energies of men and women are correlated in stable matchings. From Eq. (8) we see that $`\rho (X,Y)=\rho (X)\rho (Y)\mathrm{exp}(XY)`$ can be considered as the probability density of finding a stable matching with energies $`X`$ and $`Y`$. Due to the factor $`\left(1\mathrm{exp}(\frac{N}{X})\right)^N`$ that appears in $`\rho (X)`$ the occurrence of stable matchings with $`X>N/\mathrm{log}N`$ is very unlikely, and by symmetry, the same cutoff exists for $`Y`$. This observation is in accordance with our analysis of the GS algorithm, where we found that the energy of women in the man-optimal stable matching (i.e. the worst possible for women) is $`N/\mathrm{log}N`$, and vice versa. In the region $`\mathrm{log}N<X,Y<N/\mathrm{log}N`$ the probability distribution simplifies to $`\rho (X,Y)(XY)^{N1}\mathrm{exp}(XY)`$. It depends only on the product of men’s and women’s energies and is sharply peaked around the curve $`XY=N`$ for large values of $`N`$. Mean-value and standard deviation of the product $`XY`$ are given by

$$XY=N,\sigma (XY)=\sqrt{N}$$
(14)
where $`\mathrm{}`$ denotes the average calculated with $`\rho (X,Y)`$. At the symmetric point $`X=Y=\sqrt{N}`$ the width of the distribution is $`\mathrm{O}(1)`$ which means that relative fluctuations go to zero for $`N\mathrm{}`$.

Finally, we investigate the probability distribution of the number of stable matchings using numerical simulations. While $`S`$ has a very asymmetric distribution which decays quite slowly for large values of $`S`$, the distribution of $`\mathrm{log}S`$ can be nicely fitted with a Gaussian, as shown in Fig. 3 for $`N=50`$. Remarkably, the width of the distribution, i.e. the standard deviation $`\sigma (\mathrm{log}S)`$ approaches a constant value $`\sigma _{\mathrm{}}0.48`$ for large $`N`$ as shown in Fig. 4. We have no argument to explain this universal behavior.

## V Acceptance threshold

Up to now we have assumed that every person has to get married. In the following we soften this condition by introducing an acceptance threshold $`\mathrm{\Delta }<1`$ which is the energy of individuals who remain single. In principle, $`\mathrm{\Delta }`$ could be chosen at random for all individuals, or at least be different for men and women, but for the sake of simplicity we consider only the case where $`\mathrm{\Delta }`$ is the same for all men and all women.

Let us consider a situation where $`\mathrm{\Delta }`$ is sufficiently large that singles do not occur. The GS algorithm is readily generalized to the case with acceptance threshold. In the man-oriented version the only difference is that now women – even when unmarried – refuse proposals from men whose energy is higher than $`\mathrm{\Delta }`$. This increases the average number of proposals that are necessary to arrive at a stable matching by a factor of $`1/\mathrm{\Delta }`$. The generalization of Eq. (4) to the case $`\mathrm{\Delta }1`$ reads

$$X(\mathrm{\Delta })\frac{\mathrm{log}N+\mathrm{C}}{\mathrm{\Delta }}Y(\mathrm{\Delta })\mathrm{\Delta }\frac{N}{\mathrm{log}N+\mathrm{C}}$$
(15)
Both energies coincide for

$$\mathrm{\Delta }_\mathrm{c}=\frac{\mathrm{log}N+C}{\sqrt{N}}$$
(16)
and we conclude that for $`\mathrm{\Delta }<\mathrm{\Delta }_c`$ there exist only few or one single stable matching. This fact can be used to generate approximately sex-fair stable matchings – which is difficult to achieve by other means – using the GS algorithm and setting $`\mathrm{\Delta }=\mathrm{\Delta }_c`$.

In the limit when $`\mathrm{\Delta }`$ is sufficiently small it becomes favorable for a certain fraction of men and women to remain single. In this case it is hard to calculate the average number of proposals $`r`$ explicitly since the series of proposals is not only terminated when all men are married but also if no acceptable women are available. On the other hand the energies of men and women coincide in this regime since there exists only one stable matching. Let us assume that a woman has received $`r`$ proposals out of which she keeps only the best one. If none of the proposals are better than the threshold $`\mathrm{\Delta }`$ she remains single at the cost of $`\mathrm{\Delta }`$. Analyzing the probability distribution of the best offer (see Appendix B) yields the average energy of women

$$Y=\frac{N}{X}\left(1\mathrm{exp}(X\mathrm{\Delta })\right)$$
(17)
From the condition $`X=Y`$ we obtain

$$X^2=N\left(1\mathrm{exp}(X\mathrm{\Delta })\right)$$
(18)
In general this equation can only be solved numerically. The limiting cases are

$$X(\mathrm{\Delta })\{\begin{array}{ccc}N\mathrm{\Delta }\hfill & \mathrm{for}\hfill & \mathrm{\Delta }\frac{1}{\sqrt{N}}\hfill \\ \sqrt{N}\hfill & \mathrm{for}\hfill & \mathrm{\Delta }\frac{1}{\sqrt{N}}\hfill \end{array}$$
(19)
which means that for sufficiently small $`\mathrm{\Delta }`$ it is advantageous for almost everybody to remain single (at the energy $`\mathrm{\Delta }`$) while for $`\mathrm{\Delta }1/\sqrt{N}`$ the energies in sex-fair matchings are $`X=Y=\sqrt{N}`$, in agreement with our previous considerations.

In Fig. 5 we show the energy of men and women in the man-optimal stable matching as function of $`\mathrm{\Delta }`$ for $`N=200`$. Each data point represents the average over 1000 executions of the GS algorithm. The agreement between numerical and analytical results is very good with the exception of the cross-over region around $`\mathrm{\Delta }_c`$ where neither of the conditions that were made to treat the problem (absence of singles in the case $`\mathrm{\Delta }>\mathrm{\Delta }_\mathrm{c}`$, existence of only one stable matching in the case $`\mathrm{\Delta }<\mathrm{\Delta }_\mathrm{c}`$) is strictly fulfilled.

In the following we investigate how the number of singles depends on the acceptance threshold $`\mathrm{\Delta }`$. The probability that a woman who receives $`r`$ proposals rejects all of them and remains single is $`p_s=(1\mathrm{\Delta })^r`$. Taking into account that the number of proposals addressed to a particular woman is distributed according to a binomial distribution the average number of singles $`N_s`$ is related to the energy $`X`$ via

$$N_s(\mathrm{\Delta })=N\mathrm{exp}(X\mathrm{\Delta })$$
(20)
(see Appendix C) and can be calculated explicitly from the solution of Eq.(18). The limiting cases are

$$N_s(\mathrm{\Delta })\{\begin{array}{ccc}N\mathrm{exp}(\mathrm{\Delta }^2N)\hfill & \mathrm{for}\hfill & \mathrm{\Delta }\frac{1}{\sqrt{N}}\hfill \\ N\mathrm{exp}(\mathrm{\Delta }\sqrt{N})\hfill & \mathrm{for}\hfill & \mathrm{\Delta }\frac{1}{\sqrt{N}}\hfill \end{array}$$
(21)
When $`\mathrm{\Delta }`$ is decreased, the first single appears when $`N_s(\mathrm{\Delta })=1`$. This is the case for $`\mathrm{\Delta }=\mathrm{log}N/\sqrt{N}\mathrm{\Delta }_\mathrm{c}`$ in accordance with our assumption that there are no singles for $`\mathrm{\Delta }>\mathrm{\Delta }_c`$.

Finally, we generalize our results concerning the average number of stable matchings to the case with threshold $`\mathrm{\Delta }<1`$. The only difference compared to the calculation for $`\mathrm{\Delta }=1`$ is that the probability distribution of $`X`$ (see Appendix A) has to be replaced by

$$\rho (X)\frac{X^{N1}}{\mathrm{\Gamma }(N)}\left(1\mathrm{exp}(\frac{N\mathrm{\Delta }}{X})\right)^N$$
(22)
which takes into account that all energies $`x_i`$ have to be smaller than $`\mathrm{\Delta }`$. Making the same approximations as before we obtain for the average number of stable matchings

$$S(\mathrm{\Delta })\frac{N}{\mathrm{e}}\left(\mathrm{log}N2\mathrm{log}(\mathrm{log}N)+2\mathrm{log}\mathrm{\Delta }\right)$$
(23)
$`S(\mathrm{\Delta })`$ goes to zero at $`\mathrm{\Delta }=\mathrm{log}N/\sqrt{N}\mathrm{\Delta }_\mathrm{c}`$ in agreement with the fact that man and woman-optimal stable matchings coincide for $`\mathrm{\Delta }<\mathrm{\Delta }_\mathrm{c}`$.

## VI Unequal number of men and women

Let us now consider the more general case where the number of men and women is not the same. To be specific, we assume that there are $`N+1`$ men and $`N`$ women. It is obvious that in this case at least one man has to remain single. As before, the GS algorithm can be used to find the man-optimal and the woman-optimal stable matching, respectively, depending on who proposes and who judges. There is now a fundamental asymmetry between the man and the woman-oriented version: when women propose the algorithm stops as usual when the last woman has found a husband; in the man-oriented version, however, the series of proposals stops only when the first man has reached the bottom of his list, i.e. when he has been rejected by all women and must stay single. Therefore, the members of the minority group (the women in our case) are substantially favored.

This is illustrated in Fig. 6 where the energy of men (circles) and women (triangles) is plotted versus the acceptance threshold $`\mathrm{\Delta }`$. The data points represent averages of 1000 samples using the GS algorithm for 201 men and 200 women. Best (open symbols) and worst solutions (filled symbols) are much closer in energy now than in the symmetric case with equal number of men and women (see Fig. 5 for comparison) indicating that the total number of stable matchings is reduced. The energies in all other stable matchings lie in the two narrow corridors between the best and the worst solution. The dashed curves are $`X(\mathrm{\Delta })=X(1)\mathrm{\Delta }`$ and $`Y(\mathrm{\Delta })=Y(1)/\mathrm{\Delta }`$ obtained from the generalization of Eq. (15) to the present case while the solid curve is the solution of Eq. (18) valid for $`\mathrm{\Delta }<\mathrm{\Delta }_\mathrm{c}`$.

## VII Conclusions

We have investigated the statistical properties of stable marriages using both numerical and analytical tools. A detailed analysis of the GS algorithm is provided. This algorithm yields two extreme stable matchings and is thus particularly useful to characterize the whole set of stable matchings. The average number of stable matchings is calculated from an integral representation that can be viewed as the partition sum of two species of particles living on the unit interval, with multiplicative interaction of their center of mass coordinates. Our analytical result agrees very well with numerical simulations, in contrast to a previously obtained, less accurate formula. In addition, the fundamental relation $`XY=N`$ between the energy of men and women in stable matchings follows quite naturally from our approach. Numerical simulations indicate that the probability distribution of the logarithm of the number of stable matchings is Gaussian with a universal width $`\sigma _{\mathrm{}}0.48`$ for $`N\mathrm{}`$.

When an acceptance threshold $`\mathrm{\Delta }`$ is introduced, the number of stable matchings shrinks and singles occur. Three different regimes can be identified: For $`\mathrm{\Delta }>\mathrm{\Delta }_\mathrm{c}\mathrm{log}N/\sqrt{N}`$ there are many stable matchings whose energies extend over a wide range and no singles occur. For $`1/\sqrt{N}<\mathrm{\Delta }<\mathrm{\Delta }_\mathrm{c}`$ there are very few or only one stable matching. In this regime the energy of men and women is nearly independent of $`\mathrm{\Delta }`$ although the number of singles increases exponentially with decreasing $`\mathrm{\Delta }`$. For $`\mathrm{\Delta }<1/\sqrt{N}`$ there are many singles and only a small fraction of married individuals. When the number of men and women is not the same, the number of stable matchings is also substantially reduced. Furthermore, the members of the minority group are strongly favored compared to the majority group, even if there is only one additional member.

The general ideas discussed in this paper could become important for investigations on the labour market since nowadays information about forthcoming job opportunities, on the one hand, and qualifications of applicants, on the other hand, can be made available to everybody via the internet. Preliminary results on applications of the stable marriage problem in the context of the labour market can be found in . Further work on this subject is in progress.

## ACKNOWLEDGMENTS

We would like to thank Alain Debecker for many interesting and inspiring discussions.

## A Calculation of $`\rho (X)`$

In this appendix we calculate the distribution function $`\rho (X)`$ of the variable $`X=_ix_i`$ which is formally given by

$$\rho (X)=_0^1\mathrm{d}^Nx\delta (X\underset{i}{}x_i)$$
(A1)
According to the central limit theorem $`\rho (X)`$ converges to a normal distribution

$$\rho (X)=\frac{1}{\sqrt{2\pi \sigma ^2}}\mathrm{exp}(\frac{(X\overline{X})^2}{2\sigma ^2})$$
(A2)
with $`\overline{X}=N/2`$ and $`\sigma ^2=N/12`$ for large $`N`$. This is however only true for the central part of the distribution whereas the tails of $`\rho (X)`$ look very different. Let us introduce new variables $`q_i`$ such that $`x_i=q_iq_{i1}`$ for $`i=1,N`$ and $`q_0=0,q_N=X`$. Disregarding the constraint $`x_i<1`$ for the moment we obtain

$$\rho (X)=\frac{1}{\mathrm{\Gamma }(N)}_0^Xdq_1\mathrm{}_0^Xdq_{N1}=\frac{X^{N1}}{\mathrm{\Gamma }(N)}$$
(A3)
where the factor $`\mathrm{\Gamma }(N)`$ occurs due to the fact that the $`q_i`$ have to be arranged in ascending order to satisfy the condition $`x_i0`$. The condition $`x_i<1`$ can now be to implemented as follows: Considering the $`q_i`$ as independent random variables uniformly distributed in the interval $`[0,X]`$ we obtain a Poisson distribution $`p(x)=\lambda \mathrm{exp}(\lambda x)`$ for the distance $`x`$ between two consecutive numbers, where $`\lambda =N/X`$. The probability that all $`x_i<1`$ is therefore given by

$$P(\mathrm{all}x_i<1)=\left(_0^1dxp(x)\right)^N=\left(1\mathrm{exp}(\frac{N}{X})\right)^N$$
(A4)
neglecting the weak correlations among the $`x_i`$ due to the fact that their sum is fixed. Combining Eqs. (A3) and (A4) we obtain the final result

$$\rho (X)=\frac{X^{N1}}{\mathrm{\Gamma }(N)}\left(1\mathrm{exp}(\frac{N}{X})\right)^N$$
(A5)
valid for $`X<N/\mathrm{log}N`$.

## B Average energy of women

Let us assume that men propose in the GS algorithm and that a woman receives $`r`$ proposals out of which she keeps only the best one. If none of the offers is better than $`\mathrm{\Delta }`$ she remains single at the energy $`\mathrm{\Delta }`$. The average energy thus obtained is

$$y=\frac{1}{r+1}\left(1(1\mathrm{\Delta })^{r+1}\right)$$
(B1)
We have to take into account that the number $`r`$ of proposals that a woman receives is not fixed but distributed according to the binomial distribution

$$b(r)=\left(\genfrac{}{}{0pt}{}{R}{r}\right)p^r(1p)^{Rr}$$
(B2)
where $`R`$ is the total number of proposals that men make during the execution of the GS algorithm and $`p=1/N`$ is the probability that such an offer is addressed to a given woman. Averaging Eq. (B1) with respect to $`b(r)`$ yields

$`{\displaystyle \frac{1}{r+1}}`$ $`=`$ $`{\displaystyle \underset{r=0}{\overset{R}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{R}{r}}\right)p^r(1p)^{Rr}{\displaystyle \frac{1}{r+1}}`$ (B3)
$`=`$ $`{\displaystyle \frac{N}{R+1}}\left(1(1{\displaystyle \frac{1}{N}})^{R+1}\right)`$ (B4)
and similarly

$$\frac{(1\mathrm{\Delta })^{r+1}}{r+1}=\frac{N}{R+1}\left((1\frac{\mathrm{\Delta }}{N})^{R+1}(1\frac{1}{N})^{R+1}\right)$$
(B5)
Combining Eqs. (B1), (B4) and (B5) we obtain

$$y=\frac{N}{R+1}\left(1(1\frac{\mathrm{\Delta }}{N})^{R+1}\right)\frac{N}{R}\left(1\mathrm{exp}(\frac{R}{N}\mathrm{\Delta })\right)$$
(B6)
The average energy of women $`Y=Ny`$ is therefore related to the energy of men $`X=R/N`$ via

$$Y=\frac{N}{X}\left(1\mathrm{exp}(X\mathrm{\Delta })\right)$$
(B7)
## C Number of singles

The probability for a woman to remain single after receiving $`r`$ proposals is $`p_s=(1\mathrm{\Delta })^r`$. Averaging with respect to the binomial distribution (B2) yields

$$p_s=\underset{r=0}{\overset{R}{}}\left(\genfrac{}{}{0pt}{}{R}{r}\right)p^r(1p)^{Rr}(1\mathrm{\Delta })^r=\left(1\frac{\mathrm{\Delta }}{N}\right)^R$$
(C1)
where $`p=1/N`$ and $`R`$ is the total number of proposals made. Using $`X=R/N`$ we obtain

$$p_s=\mathrm{exp}(X\mathrm{\Delta })$$
(C2)


text:
# CLASS B0827+525: ‘Dark lens’ or binary radio-loud quasar?

## 1 Introduction

There are several indirect ways to detect the presence of dark matter in galaxies at relatively low redshifts (e.g. microlensing, rotation curves of spiral galaxies, polar-ring galaxies, etc.). However, at intermediate and high redshifts, weak and strong gravitational lensing is the only method of detecting the presence of dark matter on galaxy scales.

Gravitational lensing does not depend on the luminosity or color of the lensing mass distribution. One can therefore expect to find ‘dark-lens’ galaxies in large gravitational-lens surveys, if they make up a significant fraction of massive galaxies (e.g. Hawkins 1997). On the basis of confirmed and candidate large-separation ($``$3 arcsec) gravitational-lens systems, which were all optically selected, Hawkins (1997) concludes that around 75% of all galaxies could be ‘dark’ (i.e. extremely underluminous). This high fraction was shown to be inconsistent with observational results (Jackson et al. 1998) from the Cosmic Lens All-Sky Survey (CLASS; Browne et al. 1999; Myers et al. 1999) and also very unlikely on the basis of a number of statistical arguments (e.g. Kochanek, Falco & Muñoz 1999; Peng et al. 1999). A very high fraction of dark galaxies therefore seems to be excluded, although a much smaller fraction ($``$10%) cannot yet be ruled out. The detection of only a single dark-lens galaxy would already have severe consequences for the standard picture of galaxy formation (e.g. White & Rees 1978). Such systems would have allowed dark matter to collapse in structures massive enough to create multiple images through gravitational lensing, but at the same time prohibit baryonic matter to settle in the dark-matter halo and initiate star formation. Clearly, one has to be careful in dismissing any candidate dark-lens galaxy as being either a binary quasar or a chance alignment of unrelated quasars or AGNs, although undoubtedly for the majority of candidates this will be the case.

In this paper, we present the strongest ‘dark-lens’ candidate found in the CLASS survey: CLASS B0827+525111Not to be confused with the gravitational-lens system APM 08279+5255 (Irwin et al. 1998).. In Section 2, we present radio, optical, near-infrared and spectroscopic observations of this system. In Section 3, we derive constraints on the lens-galaxy mass-to-light ratio and compare those with observations of ‘luminous-lens’ galaxies. In Section 4, we summarize our results and suggest future work on this system.

## 2 Observations

The CLASS survey aims to find all multiply-imaged flat-spectrum ($`\alpha `$$``$0.5 for $`S_\nu `$$``$$`\nu ^\alpha `$) radio sources in the northern hemisphere with a total flux density of $`S_{5\mathrm{G}\mathrm{H}\mathrm{z}}`$$``$30 mJy, a flux-density ratio between lens images $``$10, a Galactic latitude $`|b|`$$`>`$$`10^{}`$ and a component separation $``$0.3 arcsec. The scientific goal of the survey is to create a sample of GL systems, that can be used to study the structure and evolution of lens galaxies at intermediate redshifts and constrain the cosmological parameters, in particular the Hubble parameter (H0). At present CLASS has discovered at least 17 new GL systems (Browne et al. 1999; Myers et al. 1999).

### 2.1 Radio observations

If not specified otherwise, the radio data presented in this section are flux-calibrated in the NRAO data-reduction package AIPS and self-calibrated, imaged and model-fitted in DIFMAP (Shepherd 1997). To obtain flux-densities of the two radio components of B0827+525, the uv-data was fitted by typically two Gaussian components.

As part of the CLASS sample of around 15 200 flat-spectrum radio sources observed during the 1994, 1995 and 1998 Very-Large-Array (VLA) A–array seasons (Myers et al. 1999), a 30–sec snapshot of B0827+525 was made at 8.5 GHz on 1995 August 13. The calibrated VLA image (Fig.1) shows two unresolved components that are separated by 2.8 arcsec and have a flux-density ratio of 2.7 \[i.e. $`S_\mathrm{A}`$$``$24 mJy and $`S_\mathrm{B}`$$``$9 mJy\].

With an average source density in the CLASS survey of $``$1 per sq. deg, the probability of chance alignment of two unrelated compact radio-loud quasars within $``$2.8 arcsec is around 2$`10^6`$, whereas the lensing rate is around $`10^3`$. In a sample of 15 200 sources the occurrence of such a close alignment has a probability of only 3%. B0827+525 was therefore regarded as a candidate GL system.

To improve the resolution and obtain a two-point spectral index between the radio components, a snapshot image of B0827+525 was made on 1997 January 3 with the Multi-Element Radio Linked Interferometer Network (MERLIN) at 5 GHz. A flux-density ratio of 2.5 was found \[i.e. SA$``$34 mJy and SB$``$14 mJy\] and a component separation of 2.8 arcsec, very similar to the VLA 8.5–GHz observations. Subsequent long-track MERLIN 5–GHz observations were made on 1997 December 8, with a total integration time on the source of $``$10 hours. The calibrated image (Fig.1) has an rms-noise level of 60 $`\mu `$Jy beam-1 and a resolution of $``$50 mas. Both components remain unresolved. The flux-densities are $`S_\mathrm{A}`$$``$34.0 mJy and $`S_\mathrm{B}`$$``$14.8 mJy, respectively. No sign of extended emission is detected in the image above a level of 0.6%. The similarity in flux-density ratio at 8.5 and 5 GHz strengthened the case for a lensing explanation of this system.

To improve the resolution by another order of magnitude, B0827+525 was observed with the Very-Long-Baseline-Array (VLBA) at 5 GHz on 1997 August 2. Snapshots were obtained over a range of hour angles to improve the uv-coverage. Phase referencing was done by rapid switching between the nearby strong calibrator B0828+493 (2 min) and B0827+525 (5–7 min). The total integration time on B0827+525 was 35 min. The final calibrated image (Fig.1) has a resolution of 2.5 mas and an rms-noise level of 0.1 mJy beam-1. The two radio components were each modeled by two Gaussian sub-components. Component A shows a NE–SW extension, whereas component B shows an area of weak extended emission in the NW–direction. The low SNRs of the NE–component in image A and the NW–extension of image B does not warrant modeling them by more than two subcomponents. The positions of the best-fit Gaussians and their respective flux densities are listed in Table 1.

Both radio components show fairly flat spectra from 4.9 to 8.5 GHz. The integrated flux density of B0827+525 (A+B) in the WENSS survey, (70$`\pm `$5 mJy at 325 MHz, mean epoch 1992.2), is about equal to the flux density at 1.4 GHz in the NVSS, FIRST and WSRT observations taken in the summer of 1997. This suggests that the individual component spectra must flatten to a spectral index of zero below 1.4 GHz. Alternatively, one of the component spectra could become inverted leaving room for the other to remain straight. High resolution low frequency data are needed to distinguish between these possibilities.

A 2–month monitoring campaign (8 epochs) with the WSRT at 1.4 GHz (with 15–arcsec resolution) in the period July-September 1997 indicates variability in the total flux density: the source flux density slowly decreased from 71 to 58 mJy with a typical error of about 1 mJy in individual measurements. It is not clear to which component these variations should be attributed. The variation time scale is rather short even though both components are compact. It is perhaps more likely that the 1.4–GHz variations are due to refractive interstellar scattering (RISS; e.g. Rickett et al. 1984), an explanation requiring a less extreme size for the components.

To examine the radio spectra of both components in more detail, near-simultaneous VLA 1.4, 4.9, 8.5 and 14.9–GHz (L, C, X and U–band) observations in A–array were done on 1999 June 27 (Table 2). The radio spectra for both components are shown in Fig.2. Components A and B have fairly different radio spectra with spectral indices $`\alpha _{1.4}^{8.5}=0.27\pm 0.01`$ (A) and $`\alpha _{1.4}^{8.5}=0.52\pm 0.01`$ (B) (where $`S_\nu `$$``$$`\nu ^\alpha `$). Also variability in component A at 14.9 GHz was detected from the 1998 May 15 and 1999 June 27 VLA observations, which is more likely to be intrinsic and not RISS, because the rms variability due to RISS decreases rapidly towards higher frequencies (e.g. Rickett et al. 1995; Walker et al. 1998; see also Koopmans & de Bruyn 2000).

Variability at both lower (1.4 GHz) and higher (15 GHz) frequencies makes a direct comparison of the spectral indices difficult. The formal difference in spectral index does not disqualify the source as possible GL candidate, but complicates the argumentation in a ‘dark lens’ hypothesis. We will return to this issue in Sect.4.

### 2.2 Optical and near-infrared observations

Hubble Space Telescope (HST) observations were made on 1998 April 18 with the Near-Infrared Camera and Multi-Object Spectrometer (NICMOS) at 1.6 $`\mu `$m (F160W), roughly corresponding to ground-based H–band. The NIC1 camera was used which provides a detector scale of 43 mas pixel-1 and a field-of-view of $`11^{\prime \prime }`$$`\times `$$`11^{\prime \prime }`$. The total exposure time was 2624 sec. The data were subjected to the standard NICMOS calibration pipeline in IRAF111IRAF (Image Reduction and Analysis Facility) is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy under cooperative agreement with the National Science Foundation. which (i) corrects for known instrumental effects (i.e. flat field and dark current), (ii) performs gain and flux calibration and (iii) removes cosmic rays and flags degraded and suspect data values222A detailed description of the NICMOS calibration pipeline process, zero-point fluxes and PSF can be found in the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) handbook for Cycle 10, available at the STScI web page.. The image was cut around B0827+525 and rotated to correct to standard orientation. The final calibrated image is shown in Fig.3 (middle). The primary northern component (A) is unresolved and clearly detected ($`m_{\mathrm{F160W}}`$=19.6$`\pm `$0.2 mag). Component B is just barely detected as a region of extended emission ($`m_{\mathrm{F160W}}`$=22.6$`\pm `$0.2 mag). No sign of a lensing galaxy is seen in the image to $`m_{\mathrm{F160W}}`$=23 mag, although the extended emission near image B could partly be from a possible lens galaxy. No optical correction (i.e. deconvolution) was applied to the image, because the PSF distortion222A detailed description of the NICMOS calibration pipeline process, zero-point fluxes and PSF can be found in the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) handbook for Cycle 10, available at the STScI web page. is not expected to change appreciably over the image separation of only 2.8 arcsec.

Similarly, B0827+525 was observed in K–band on 1998 May 9 with the Near InfraRed Camera (NIRC) at Keck. A total of 25 frames, each having a 60–sec integration time, were obtained of the field. The telescope was moved by approximately 10–15 arcsec between individual integrations. The seeing at the time of the observations was 0.6–0.75 arcsec (FWHM). The exposures were combined and rotated to correct the orientation and the resulting image is shown in Fig.3 (right). The observations confirm the conclusions drawn from the HST exposure, that image A is consistent with a point source and image B seems somewhat extended.

An R–band image was taken on 1998 May 19 on the WHT. Two 600-sec exposures were obtained with a seeing of 1.4–1.5 arcsec (FWHM). It was re-observed on 1998 May 20 with two exposures of 600–sec and one of 300–sec exposure, respectively. The seeing was around 1.0 arcsec (FWHM). All exposures were combined and reduced in the standard way in the data-reduction package IRAF. The resulting image (Fig.3, left) only show image A, whereas image B is not detected above the noise level. Unfortunately no useful calibration was obtained, but based on the exposure time and expected noise-level, we roughly estimate for component A an apparent magnitude of $`m_\mathrm{R}`$$``$22.0–22.5 mag and for component B a lower limit $`m_\mathrm{R}`$$``$24.0 mag.

### 2.3 Spectroscopic observations

B0827+525 was observed with the Low Resolution Imaging Spectrograph (LRIS; Oke et al. 1995) on the W. M. Keck II Telescope on the night of 1998 April 21. Three exposures were taken in longslit mode with a total exposure time of 4500 sec. The 300 grooves mm-1 grating was used, giving a dispersion of 2.44 Å per pixel. The slit was placed along both components, with a position angle of 50. The final spectrum covers a wavelength range of 4024–9012 Å. The spectra were reduced using standard IRAF tasks. A correction for the response of the CCD was determined from observations of the Oke standard star BD332642 (Oke 1990). The corrected spectra were weighted by the squares of their SNRs and co-added to produce the final B0827+525 spectrum. Fig.4 shows the spectrum of component A, which has been smoothed using a boxcar kernel with a width of five pixels. Emission lines associated with Si-IV $`\lambda 1397`$, C-IV $`\lambda 1549`$, He-II $`\lambda 1641`$, and possibly C-III\] $`\lambda 1909`$ are seen, establishing a source redshift of 2.064. No emission at the position of component B is detected. We also find no evidence in the spectrum of component A for a second redshift, possibly from a lens galaxy.

## 3 Analysis

Before proceeding with our analysis, let us summarize the most important observational results. We have found that B0827+525 consists of two radio components, separated by 2.8 arcsec. At least component A is variable, and both components have substructure at a scale of a few mas. The integrated radio spectra of both components from 1.4 to 14.9 GHz are somewhat different. The redshift for the brightest radio component (A) is 2.064. Optical counterparts for both radio components are found in H and K–band observations. The emission associated with the fainter radio component appears extended, which could suggest that it is very faint galaxy emission.

### 3.1 The ‘dark lens’ hypothesis

To compare the observed limits on the mass-to-light ratio of a possible lens galaxy with those found for typical lens galaxies, we first calculate the velocity dispersion of the lensing mass, using the Singular-Isothermal-Sphere (SIS) mass model (Binney & Tremaine 1987). This model relates the lens-image separation ($`\mathrm{\Delta }\theta `$) to a line-of-sight velocity dispersion ($`\sigma `$), given the cosmological model and the source- and lens redshifts:

$$\mathrm{\Delta }\theta =8\pi \left(\frac{\sigma }{c}\right)^2\left(\frac{D_{\mathrm{ds}}}{D_\mathrm{s}}\right),$$
(1)
where $`c`$ is the velocity of light, and $`D_{\mathrm{ds}}`$ and $`D_\mathrm{s}`$ are the angular-diameter distances between lens-source and observer-source, respectively (e.g. Schneider, Ehlers and Falco 1992). The angular-diameter distance is a function of redshifts and cosmological model. Using the component separation, we can furthermore calculate the mass contained inside the Einstein radius

$$M(\mathrm{\Delta }\theta )=3.110^{10}\times \left(\frac{\mathrm{\Delta }\theta }{1^{\prime \prime }}\right)^2\left(\frac{D}{\mathrm{Gpc}}\right)\text{ M}_{},$$
(2)
where $`D`$$``$$`D_\mathrm{d}D_\mathrm{s}/D_{\mathrm{ds}}`$ and $`D_\mathrm{d}`$ is the angular-diameter distance to the lens galaxy. To calculate the luminosity of the lens galaxy, using the filter band $`\lambda `$ (e.g. U, B, V, $`\mathrm{}`$), we subsequently use the relation

$$L_\lambda /\mathrm{L}_{\lambda ,}=10^{0.4(M_{\lambda ,}m_\lambda +\mathrm{DM}+\mathrm{K}_\lambda )},$$
(3)
where DM is the distance modulus of the lens galaxy, Mλ,⊙ is the absolute magnitude of the sun, $`m_\lambda `$ is the apparent magnitude of the lens galaxy and Kλ is the K–correction for the filter band $`\lambda `$. We use the lower limit of $`m_{\mathrm{F160W}}`$$`>`$23 mag on the lens galaxy, found from the HST exposure (Fig.3) to constrain the mass-to-light ratio in H–band. We use the H–band observations, because K corrections are relatively small ($`<`$0.2 mag for $`z`$$``$1 and $`<`$0.5 mag for $`z`$$``$2; Poggianti 1997). We do not use evolutionary corrections and assume $`M_{\mathrm{H},}`$=3.46 mag.

For a source redshift of 2.064, we have plotted the H–band mass-to-light ratio in Fig.5 for two different cosmologies (A: $`\mathrm{\Omega }_\mathrm{m}`$=1 and $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0; B: $`\mathrm{\Omega }_\mathrm{m}`$=0.3 and $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7) and H0=65 km s-1 Mpc-1. We assume a Friedmann-Robertson-Walker universe. A minimum mass-to-light ratio around 100 is found, using cosmology B. This mass-to-light ratio implies a lens-galaxy redshift around 1.5 and a high velocity dispersion of 400–500 km/s. For cosmology A, the minimum mass-to-light ratio would increase to about 200. From Jackson et al. (1998) we find that the mass-to-light ratio in H–band for B0827+525 is a factor 15–100 larger than the H–band mass-to-light ratios of typical lens galaxies in the CLASS survey.

In Fig.6, we have plotted the expected H–band magnitude of the lens galaxy as function of redshift, galaxy type and cosmology. We used (i) the velocity dispersion determined from eqn.1, (ii) the relation between velocity dispersion and B–band magnitude of elliptical and spiral galaxies from Fukugita et al. (1991) and (iii) the colors and K–corrections from Poggianti (1997). From Fig.6 an upper limit on $`m_\mathrm{H}`$ of 17 mag is found, nearly independent of galaxy type. If evolutionary corrections are applied (Poggianti 1997), we find no significant differences in these results. The upper limit is at least 6 mag brighter than the lower limit on the H–band magnitude of the possible lens galaxy, whereas for all confirmed GL systems (Jackson et al. 1998), for which the lens and source redshifts are known, differences between the observations and the model are $``$1 mag.

Hence, if B0827+525 is a lens system then the object is indeed a ‘dark lens’ system. However, also image B must be darkened and reddened by the ‘dark lens’ galaxy. For further discussion of this possibility we refer to Sect.4.

### 3.2 The binary quasar hypothesis

As argued in Sect.2.1 the probability of a chance alignment of two radio-loud quasars is rather small. Hence, if the source is not a lens systems it is most likely a pair of physically close quasars. We call this the ‘binary quasar’ hypothesis even though we have no evidence for the quasar nature of component B in the optical.

The compact flat-spectrum radio source population selected via CLASS is generally believed to emit its radio emission in a highly non-isotropic manner (e.g. Orr & Browne 1982) with the radio axis pointed in our direction within a small angle. As first argued by Scheuer & Readhead (1979) there should then be many radio-quiet quasars (QSO) for every radio-loud quasar (QSR). It is then obvious that the probability of finding two related radio-loud quasars pointed at us is rather small. In fact we can probably state that for every QSR–QSR pair there should be around 20 QSR-QSO pairs in the CLASS survey, if $``$5% of all QSOs are also radio-loud, i.e. $``$30 mJy at 8.5 GHz (e.g. Hooper et al. 1995; Bischof & Becker 1997).

Under the assumption that B0827+525 is the only QSR–QSR pair (i.e. binary radio-loud quasar) in CLASS, the probability that a flat-spectrum radio source is part of a QSR–QSR pair is around 1 in 15 000. This means that the probability that a CLASS source is part of a QSR–QSO pair is around 1 in 750. In the Large Bright Quasar Survey (LBQS; Hewett et al. 1998) two QSO–QSO pairs were found from a sample of around 1000 optically selected quasars, hence 1 in 500. This number is very close to that found from the CLASS survey, which means that the presence of one QSR–QSR pair in the CLASS survey is consistent with the rate of QSO–QSO pairs in optical quasar surveys (see also Kochanek et al. 1999).

However, if B0827+525 is a binary radio-loud quasar and we compare it with the list of wide separation quasar pairs in Kochanek et al. (1999), we notice two things: First, B0827+525, has the smallest separation (2.8 arcsec or 23 kpc/$`h_{50}`$) of all quasar pairs. Second, only 2 out of 13 non-lens quasar pairs have higher redshifts (i.e. LBQS 1429-008 and Q2345+007). On the other hand, B0827+525 would also be one of the largest separation lens systems (e.g. Browne et al. 1999; Myers et al. 1999; CASTLES survey222http://cfa-www.harvard.edu/glensdata/). B0827+525 thus appears to be at the parameter-space border delineated by the optical quasar pairs listed by Kochanek et al. (1999) and the confirmed GL systems. It is difficult to see how significant these issues are in the context of binary quasars and how it possibly relates to the fact that both quasars are radio-loud.

## 4 Discussion & Conclusions

If taken at face value, most of the observational evidence seems to be in favor of B0827+525 being a binary radio-loud quasar: (i) somewhat different radio spectra from 1.4 to 14.9 GHz for the two components, (ii) a different radio and optical brightness ratio, (iii) the weaker radio component appears slightly resolved, (iv) the extended nature of the optical and near-infrared emission near component B, whereas that of image A is compact and (v) the agreement with the number statistics of binary radio-loud quasars in the CLASS survey with those of optical surveys.

However, the most interesting explanation of B0827+525, a “dark lens” system, cannot convincingly be excluded and in view of the cosmological importance of the presence of galaxy-sized concentrations of dark matter – even though they must be rare (e.g. Jackson et al. 1998) – we now discuss why we believe the ‘dark lens’ hypothesis can not yet be discarded: (1) We have seen that at least one of the radio components is variable. If they are lens images, a time delay between them combined with variability of the source can result in a difference in simultaneously measured spectral indices, especially if the time delay is of the order of the variability time scale. (2) If we place a massive ‘dark lens’ galaxy near component B, containing a large amount of dust, this galaxy would obscure most optical and a significant fraction of the near-infrared emission from component B. In the secure lens system B0218+357 a similarly large discrepancy has been observed between the optical and radio brightness ratios of the lens images (e.g. CASTLES survey), which can also be explained in terms of obscuration by the high-column density of the rich ISM in the lens galaxy (e.g. Wiklind & Combes 1995; Menten & Reid 1996; Combes Wiklind & Nakai 1997; Combes & Wiklind 1997; Gerin et al. 1997; Combes & Wiklind 1998). This dust, however, would probably not block all of the lens-galaxy emission. In fact, most of the extended emission seen near radio component B could be coming from the lens galaxy itself and not from the quasar. (3) The optical brightness of image B is surprisingly low for a radio-loud quasar, which might suggest some form of extinction. (4) A high scattering measure, associated with the ionized component of the ISM in the lens galaxy (see for example Marlow et al. 1999; Jones et al. 1996), could furthermore scatter-broaden radio component B, explaining why it appears somewhat resolved.

So far in the ‘dark lens’ hypothesis we have assumed that the lensing mass is a single galaxy dominated by non-baryonic matter in its inner parts and therefore extremely underluminous. If the lensing mass distribution is a high-redshift ($`z`$$``$1) group or cluster, the constraints on the mass-to-light ratio, which is usually an order of magnitude larger in clusters and groups of galaxies than in the inner parts of galaxies (e.g. van der Marel 1991; Carlberg, Yee & Ellingson 1997), is somewhat alleviated. However, for most of the lens galaxies from the CLASS survey, the H–band mass-to-light ratios are around unity (Jackson et al. 1998). This is two orders of magnitude smaller than the minimum H–band mass-to-light ratio found in Sect.3.1 and can therefore not entirely account for the high mass-to-light ratio within the context of a ‘normal’ galaxy evolution. Hence also in the case that the lensing mass is a group or cluster, it must be underluminous. A possible candidate for this type of mass concentration was recently found by Erben et al. (2000).

Thus, although the binary quasar hypothesis seems more likely at face value, the ‘dark lens’ hypothesis can not be ruled out convincingly. However, the latter hypothesis does give concrete predictions: First, the substructure in both radio components should be related, although distorted by the intervening gravitational potential. We have obtained multi-frequency global VLBI data on B0827+525 at three frequencies to test this. (2) If the variability of the radio components is intrinsic, one should be able to correlate their light-curves to find a time-delay. This could definitively prove or disprove whether B0827+525 is a ‘dark lens’ system. (3) If the extended H and K–band emission seen near component B is partly from a lens galaxy, it should have a redshift smaller than 2.064. Near-infrared spectroscopy of this emission should therefore be attempted.

## Acknowledgments

The authors would like to thank Lee Armus and David Hogg for obtaining a Keck NIRC K–band image. LVEK and AGdeB acknowledge the support from an NWO program subsidy (grant number 781-76-101). This research was supported in part by the European Commission, TMR Programme, Research Network Contract ERBFMRXCT96-0034 ‘CERES’. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Westerbork Synthesis Radio Telescope (WSRT) is operated by the Netherlands Foundation for Research in Astronomy (ASTRON) with the financial support from the Netherlands Organization for Scientific Research (NWO). MERLIN is a national UK facility operated by the University of Manchester on behalf of PPARC. This research used observations with the Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by Associated Universities for Research in Astronomy Inc. under NASA contract NAS5-26555.


text:
# NONPERTURBATIVE CALCULATION OF SCATTERING AMPLITUDES11footnote 1To appear in the proceedings of the fourth workshop on Continuous Advances in QCD, Minneapolis, Minnesota, May 12-14, 2000.

## 1 Introduction

Recently there has been a significant amount of progress in the nonperturbative calculation of bound-state properties in light-cone quantized field theories.  Comparable progress in the calculation of scattering states has been lacking, although techniques have been proposed.  The ability to compute outcomes of scattering processes is not only of interest in its own right, but also useful for construction of renormalization conditions in bound-state calculations.  Here we consider a calculational scheme based on numerical approximations to resolvents  and on a generalization of a Hamiltonian formulation for the scattering of composites.

The construction is done in terms of light-cone coordinates, where $`t+z`$ takes the role of time.  This is not a necessary choice, but it does bring several advantages and makes immediate contact with work that uses these coordinates for bound-state problems.  The advantages include a simple vacuum, well defined Fock-state expansions, and kinematical boosts. This last advantage is particularly important for the scattering of composites because one can compute the state of a composite particle in one frame and boost to other frames without again referencing the interaction.

The remainder of this paper is structured as follows: Section 2 contains a description of the notation used for light-cone coordinates and a formulation of two-body cross sections in terms of these coordinates. The key quantity in this formulation is a matrix element of a resolvent; a numerical method for its computation is discussed in Sec. 3 and applied to a nonrelativistic example in Sec. 4. An extension to the scattering of composites is given in Sec. 5. Finally, Sec. 6 provides a summary and indications of future work.

## 2 Cross Sections in Light-Cone Coordinates

We use the following notation for light-cone coordinates:

$$x^\pm =t\pm z,𝐱_{}=(x,y),$$
(1)
where $`x^+`$ is taken as the time variable. Momentum components are similarly defined

$$p^\pm =E\pm p_z,𝐩_{}=(p_x,p_y).$$
(2)
A dot product is then written as

$$px=\frac{1}{2}(p^+x^{}+p^{}x^+)𝐩_{}𝐱_{},$$
(3)
so that $`p^{}`$ is conjugate to $`x^+`$ and is recognized as the light-cone energy. The other momentum components form the light-cone analog of the three-momentum

$$\underset{¯}{p}(p^+,𝐩_{}).$$
(4)
The evolution of a state is determined by the Hamiltonian $`𝒫^{}`$; however, the term “light-cone Hamiltonian” is also used  for the combination

$$H_{\mathrm{LC}}=P^+𝒫^{},$$
(5)
where $`P^+`$ is the total momentum. The fundamental eigenvalue problem is

$$H_{\mathrm{LC}}|p=\frac{M^2+p_{}^2}{p^+/P^+}|p,\underset{¯}{𝒫}|p=\underset{¯}{p}|p.$$
(6)
For two-body scattering $`A+BC+D`$, the center-of-mass cross section can be written

$$\frac{d\sigma }{d\mathrm{\Omega }}_{\mathrm{cm}}=\frac{1}{2E_A2E_Bv_{\mathrm{rel}}}\frac{|\stackrel{}{p}_C||_{fi}|^2}{16\pi ^2E_{\mathrm{cm}}},$$
(7)
where $`M_{fi}`$ is the invariant amplitude extracted from the $`S`$ matrix

$$S_{fi}=f|i+(2\pi )^4\delta ^{(4)}(p_fp_i)i_{fi},$$
(8)
with states normalized according to

$$f|i=(16\pi ^3)^2p_A^+\delta (\underset{¯}{p}_A\underset{¯}{p}_C)p_B^+\delta (\underset{¯}{p}_B\underset{¯}{p}_D).$$
(9)
The flux factors are determined by the current matrix element

$$A|J^0|A=\frac{1}{2}A|(J^++J^{})|A=2E_A,$$
(10)
with $`E_A`$ the energy of particle $`A`$. For fermions, the spinors are normalized by $`\overline{u}u=2m`$. From the optical theorem we also have

$$\sigma =\frac{2m_{ii}}{2E_A2E_Bv_{\mathrm{rel}}}.$$
(11)
In the light-cone interaction picture  the invariant amplitude can be related to a matrix element of a $`T`$ operator. The $`S`$ matrix can be written as

$$S_{fi}=f|U(x^+,x_0^+)|i=f|i2\pi iT_{fi}\delta (p_f^{}p_i^{}),$$
(12)
with the $`T`$ operator defined by

$$T=V+V\frac{1}{p^{}+iϵ𝒫^{}}V,$$
(13)
where

$$𝒫^{}=𝒫_0^{}+V.$$
(14)
On removal of the momentum conserving delta function

$$T_{fi}=t_{fi}\delta (\underset{¯}{p}_f\underset{¯}{p}_i),$$
(15)
we obtain

$$_{fi}=\frac{1}{16\pi ^3}t_{fi}.$$
(16)
We then need a means to compute $`t_{fi}`$. In particular, we need to handle matrix elements of the second term in $`T`$ which involve the resolvent $`(p^{}+iϵ𝒫^{})^1`$.

## 3 Recursion Method for Resolvents

We divide the matrix element of the $`T`$ operator into two pieces, the Born term and the resolvent term

$$T_{fi}=V_{fi}+T_{fi}^{(2)},T_{fi}^{(2)}f|V\frac{1}{p^{}+iϵ𝒫^{}}V|i.$$
(17)
To the latter we apply the recursion method developed by Haydock  in the context of condensed matter physics. This method is based on the Lanczos diagonalization algorithm  which converts the original Hamiltonian operator into an approximate tridiagonal representation. The tridiagonal form is easily inverted to yield the resolvent.

Define an initial vector $`|u_1=cV|i`$, with $`c`$ a normalization constant, and carry out Lanczos recursions

$`b_1`$ $`=`$ $`0,a_n=u_n|𝒫^{}|u_n,u_n|u_n=1,`$
$`b_{n+1}|u_{n+1}`$ $`=`$ $`𝒫^{}|u_{n+1}a_n|u_nb_n|u_{n1}.`$ (18)
Although this could in principle be done analytically, the usual approach is to begin with a large but sparse matrix representation for $`𝒫^{}`$ in a Fock basis. A standard way of constructing such a representation is discrete light-cone quantization (DLCQ).  The vectors $`|u_n`$ provide an orthonormal basis in which $`𝒫^{}`$ has the tridiagonal representation

$$𝒫^{}\left(\begin{array}{cccccc}a_1\hfill & b_2\hfill & 0\hfill & 0\hfill & 0\hfill & \mathrm{}\hfill \\ b_2\hfill & a_2\hfill & b_3\hfill & 0\hfill & 0\hfill & \mathrm{}\hfill \\ 0\hfill & b_3\hfill & a_3\hfill & b_4\hfill & 0\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & b_4\hfill & .\hfill & .\hfill & \mathrm{}\hfill \\ 0\hfill & 0\hfill & 0\hfill & .\hfill & .\hfill & \mathrm{}\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & \mathrm{}\hfill \end{array}\right)A,$$
(19)
and the forward matrix element of the resolvent is approximated by the upper left matrix element of the inverse

$$T_{ii}^{(2)}\frac{1}{c^2}\left\{[(p_i^{}+iϵ)IA]^1\right\}_{1,1}.$$
(20)
Non-forward matrix elements can be obtained by using the block Lanczos procedure in place of the ordinary one described above.

Unfortunately, the discrete nature of the calculation does not allow this finite matrix representation to be useful in all respects. In the limit of $`ϵ0`$, the imaginary part of this approximation to $`T_{ii}^{(2)}`$ will contain only delta functions at the discrete eigenvalues of the finite matrix approximation $`A`$. To properly model the contribution from the continuum, an extension must be made. This can be done in terms of a continued fraction representation

$$c^2T_{ii}^{(2)}=\frac{1}{p_i^{}+iϵa_1\frac{b_2^2}{p_i^{}+iϵa_2\frac{b_3^2}{p_i^{}+iϵa_3\mathrm{}}}},$$
(21)
by replacing the last level with a terminator, a model for an infinite fraction.

A model can be generated in various ways; however, a logical choice in the present situation is to use the second Born approximation as a generator. The matrix element of the resolvent for the free Hamiltonian

$$i|V\frac{1}{p^{}+iϵ𝒫_0^{}}V|i$$
(22)
is computed exactly and then in continued fraction form by iterating with $`𝒫_0^{}`$ an equal number of times. The exact form dictates the structure of the remaining infinite fraction, which is used as a terminator for the full resolvent.

## 4 Nonrelativistic Example

As an example of how the terminated recursion method works, we consider nonrelativistic S-wave scattering from an exponential potential. This is a common choice for testing numerical methods for scattering problems because it has an exact solution.

The partial wave expansion for the $`T`$-matrix is

$$\stackrel{}{k}^{}|T(E)|\stackrel{}{k}=\frac{1}{2\pi ^2}\underset{l}{}(2l+1)T_l(k^{},k;E)P_l(\widehat{k}^{}\widehat{k}),$$
(23)
where

$`T_l(k^{},k;E)`$ $`=`$ $`{\displaystyle 𝑑ru_{k^{}l}(r)V(r)u_{kl}(r)}`$ (25)
$`+{\displaystyle 𝑑ru_{k^{}l}(r)V(r)\frac{1}{E+iϵH_l}V(r)u_{kl}(r)}`$
and

$$u_{kl}(r)=rj_l(kr),H_l=\frac{1}{2m}\frac{d^2}{dr^2}+\frac{l(l+1)}{2mr^2}+V(r).$$
(26)
The partial wave contribution $`T_l`$ is related to the phase shift $`\delta _l`$ by

$$T_l(k,k;k^2/2m)=\frac{1}{2mk}e^{i\delta _l}\mathrm{sin}\delta _l\text{or}\mathrm{tan}\delta _l=\frac{mT_l}{eT_l}.$$
(27)
For S-wave scattering of a particle of mass $`m`$ from $`V(r)=e^{r/a}/ma^2`$
the exact solution is

$$\delta _0(k)=\text{Arg}(J_{2ika}(2\sqrt{2})\mathrm{\Gamma }(1+2ka))2ka\mathrm{ln}\sqrt{2}.$$
(28)
To solve this problem numerically, 333The method presented here is not by any means the best way to solve this particular problem. It is introduced purely for the purpose of testing the recursion method. we introduce a discrete approximation with $`M`$ interior points between $`r=0`$ and a radial cutoff $`R`$. The step size is

$$\mathrm{\Delta }r=\frac{R}{M+1}.$$
(29)
Integrals are approximated by the trapezoidal rule, and derivative operators by finite differences. The radial cutoff is kept larger than $`2\pi /k`$ in order to include one complete wavelength at the chosen energy. The terminator for the recursion is generated from the second Born approximation

$$T_l^{(2)}𝑑r𝑑r^{}u_{k^{}l}(r^{})V(r^{})g_{El}^{(0)}(r^{},r)V(r)u_{kl}(r),$$
(30)
where

$$g_{El}^{(0)}(r^{},r)=\frac{mi\kappa }{2\pi }rr^{}j_l(\kappa r_<)h_l^+(\kappa r_>)$$
(31)
and $`\kappa =\sqrt{2mE}`$.

Figure 1 displays the results of calculations done for $`k=0.25/a`$, $`R=30a>2\pi /k25a`$, and $`M=49`$, 99, 199, 399, and 799. They show good convergence properties. More importantly, they show that the imaginary part of the $`T`$-matrix element is calculated reliably.

## 5 Scattering of Composites

The formulation in Sec. 2 assumes implicitly that the particles involved in the scattering process are elementary. In a field-theoretic setting this is not typically the case for a nonperturbative calculation. The particles $`A`$, $`B`$, $`C`$, and $`D`$ can be themselves complicated bound states of the full Hamiltonian. It is then difficult to determine the nature of the interaction $`V`$ that appears in the $`T`$ operator.

To construct the form of the $`T`$ operator in this case, we generalize a construction given by Wick  to compute the $`S`$ matrix

$$S_{ABCD}=CD|AB+$$
(32)
from scattering eigenstates of the light-cone Hamiltonian $`H_{\mathrm{LC}}`$. Let $`A^{}`$ and $`B^{}`$ create one-particle solutions from the Fock vacuum $`|0`$, so that

$$|A=A^{}|0,H_{\mathrm{LC}}|A=\frac{m_A^2+p_A^2}{p_A^+/P^+}|A.$$
(33)
The two-particle scattering eigenstates of $`H_{\mathrm{LC}}`$ are then

$$|AB\pm =B^{}|A+\frac{1}{s_{AB}\pm iϵH_{\mathrm{LC}}}V_B|A,$$
(34)
where

$$V_B[H_{\mathrm{LC}},B^{}]\frac{m_B^2+p_B^2}{p_B^+/P^+}B^{}$$
(35)
and

$$s_{AB}=\frac{m_A^2+p_A^2}{p_A^+/P^+}+\frac{m_B^2+p_B^2}{p_B^+/P^+}$$
(36)
is the invariant mass of the scattering state and the Mandelstam $`s`$ of the scattering process. The definition of $`V_B`$ is chosen to ensure that

$$H_{\mathrm{LC}}|AB\pm =s_{AB}|AB\pm .$$
(37)
The $`T`$ operator is constructed by reducing the $`S`$ matrix to

$$S_{ABCD}=\delta _{CD,AB}2\pi i\delta (s_{AB}s_{CD})T_{\mathrm{LC}fi}.$$
(38)
We find

$$T_{\mathrm{LC}fi}=P^+T_{fi}=CD|V_B|A,$$
(39)
which can be expanded as

$$P^+T_{fi}=C|V_D^{}\frac{1}{s_{AB}+iϵH_{\mathrm{LC}}}V_B|A+C|DV_B|A,$$
(40)
with

$`C|DV_B|A`$ $`=`$ $`C|V_B{\displaystyle \frac{1}{s_{AB}\frac{m_B^2+p_B^2}{p_B^+/P^+}\frac{m_D^2+p_D^2}{p_D^+/P^+}H_{\mathrm{LC}}}}V_D^{}|A`$ (41)
$`+C|[D,V_B]|A.`$
The first term in (40) can be interpreted as a direct scattering of $`B`$ by $`A`$. The first term in (41) is the crossed contribution, where $`D`$ is produced before $`B`$ is absorbed. The second term in (41) is a generalized seagull contribution, which did not appear in Wick’s formalism.

To use these expressions efficiently, $`C|DV_B|A`$ should be computed explicitly rather than work with the decomposition into crossed and seagull terms. This is because the decomposition introduces an inverse of the Hamiltonian. The inverse is unavoidable for the direct term, and here one would apply the recursion method discussed in the previous sections. For both terms one needs the action of $`V_B`$ on the one-particle state $`|A`$.

## 6 Concluding Remarks

A method for the nonperturbative calculation of cross sections has been presented in the context of light-cone quantized field theories. The key expression is (40), where the light-cone $`T`$ matrix is given in terms of known operators and composite-particle eigenstates. The terminated recursion technique of Haydock  is to be used to handle the inverse of the light-cone Hamiltonian. Work on an application to a simple model  is in progress. One direction to consider for possible improvements to the method is to explore variational formulations analogous to the Schwinger and Kohn variational principles  used in nonrelativistic scattering.

## Acknowledgments

The work reported here was supported in part by the Department of Energy, contract DE-FG02-98ER41087.

## References