text:
# The Quintuple Quasar: Radio and Optical Observations

## 1 Introduction

The quasar PMN J0134–0931 calls attention to itself in many ways. It is very bright at radio frequencies (1 Jy at 2 GHz); its optical counterpart is extraordinarily red ($`BK11`$); and it has a unique radio morphology consisting of six compact components within a circle of diameter $`0\stackrel{}{\mathrm{.}}7`$. It is therefore unsurprising, in retrospect, that it was discovered to be an interesting object by two independent groups.

Winn et al. (2002a; hereafter, W02) found it during a survey for radio-loud gravitational lenses. They showed that at least five of the six radio components have identical continuum radio spectra, and discovered a curved arc of radio emission joining two of the components. These properties are hallmarks of strong gravitational lensing.

Gregg et al. (2002; hereafter, G02) found it in a cross-compilation of radio, optical, and near-infrared catalogs that was designed to identify the reddest quasars on the sky. They measured its redshift ($`z_s=2.2`$) and found at least three components in a near-infrared image. The multiplicity of the near-infrared counterpart, and the extremely large inferred luminosity ($`M_R=29.6`$, even under the conservative assumption of no extinction) were evidence of gravitational lensing.

However, neither group could explain the unique morphology with a plausible lens model, because strong lensing by a single galaxy usually results in only two or four images. Another difficulty was that two of the components were found to have a lower surface brightness than the other components in high-resolution radio images, in apparent contradiction of the conservation of surface brightness by gravitational lensing.

Also unknown was the reason for the extremely red color of J0134–0931. The two most obvious possibilities—reddening due to dust in the host galaxy, or due to dust in the lens galaxy—could not be distinguished. Hall et al. (2002) detected strong Ca II H and K absorption lines at $`z_l=0.76`$ (presumably the lens redshift), but they were not able to determine the reason for the redness of the quasar spectrum.

In a companion paper, Keeton & Winn (2003; hereafter, KW03) present the first quantitative models that explain all the previously known properties of J0134–0931. In those models, five of the radio components are multiple images of a single quasar, and the sixth component represents a different source, which is presumably a second component of the same background radio source. The positions of the two lens galaxies are constrained remarkably well by the observed image configuration. Furthermore, the two low-surface-brightness components are located close to the expected position of one of the lens galaxies, suggesting that scatter broadening by ionized material in that lens galaxy could account for their large angular sizes.

In addition to explaining the previous data, the KW03 models make a number of predictions about more subtle properties of J0134–0931 that would be detectable in radio and optical observations with high sensitivity and high angular resolution. In this paper we present multi-frequency observations of J0134–0931 with the Very Long Baseline Array (VLBA111The Very Long Baseline Array (VLBA) is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities for Research in Astromomy, Inc.) and the Hubble Space Telescope (HST222Data from the nasa/esa Hubble Space Telescope (HST) were obtained from the Space Telescope Science Institute, which is operated by aura, Inc., under nasa contract NAS 5-26555.) that confirm some of these predictions, and refute none of them.

Before describing the new observations, in § 2 we provide an overview of the previously-known radio morphology of J0134–0931, and the nomenclature used to describe it. We also review the main features of the KW03 models. The radio data (§ 3) confirm that the sixth and faintest radio component has a significantly different spectral index than the other five components (§ 3.3). The two components that were observed to have larger angular sizes are indeed being scatter broadened (§ 3.4). There is tantalizing (but tentative) evidence for one of the additional images of the second background source that are predicted by the KW03 models (§ 3.5). The extended radio emission (§ 3.6) and component shapes (§ 3.7) can be explained by the models.

After subtracting point sources from the optical images (§ 4.1), we identify two faint but significant peaks of residual flux as direct detections of the lens galaxies (§ 4.2). The optical counterparts of the radio components have different colors, and their optical flux ratios differ from their radio flux ratios, implying that the extreme red color of J0134–0931 is at least partly due to reddening by dust in the lens galaxy (§ 4.3). Finally, in § 5, we summarize the evidence that J0134–0931 is the first known instance of a quintuply-imaged quasar.

## 2 Overview

Figure 1 is a diagram of the observed configuration of the radio components. The dots represent the six compact radio components A–F, which are named in decreasing order of radio flux density. The size of each dot is proportional to the logarithm of its flux density ratio relative to component D. The angular sizes of C and E are much larger than the other components at 1.7 GHz and 5.0 GHz, which is why we have used gray dots to represent those components instead of black dots. There is a curved arc extending between A and B at 1.7 GHz, represented schematically by a solid line. Components A–E all have the same spectral index, as measured from 5 GHz to 43 GHz ($`\alpha =0.69\pm 0.04`$, where $`S_\nu \nu ^\alpha `$), but the spectral index of F was unknown.

According to the KW03 models, the background radio source has two components (S1 and S2) and there are two lens galaxies (Gal-N and Gal-S). The five observed radio components A–E are five images of S1, and F is the brightest of three images of S2. The positions and masses of the lens galaxies are fairly well constrained by the data, but their ellipticities and position angles are not. For reference, the critical curves and caustics of a representative lens model are shown in Fig. 2. The reader is referred to KW03 for a full theoretical analysis.

## 3 Radio evidence

### 3.1 Observations and calibration

This system was observed with the VLBA at 1.7 GHz and 5.0 GHz by W02. Below, we present the results from new VLBA observations at three different frequencies. We also present images based on both the old and the new data, to provide a uniform multi-frequency study. A list of the basic parameters of all the VLBA observations is given in Table 1.

The new observations were conducted in standard VLBA continuum frequency bands at 2.3 GHz, 8.4 GHz, and 15 GHz. The array consisted of the usual ten antennas, except that the Pie Town antenna did not participate in the 15 GHz observation.333Instead, one antenna from the Very Large Array (VLA) participated. However, the data from this antenna were much noisier than the VLBA data, and we did not include data from the VLA antenna in our final analysis. For the 2.3 GHz observation, we did not employ phase referencing because the source was sufficiently bright for self-calibration. At the higher frequencies, we switched between J0134–0931 and the nearby, bright, compact source J0141–0928 in order to calibrate the relative antenna gains. The cycle times were 4 minutes and 2 minutes at 8.4 GHz and 15 GHz, respectively. Short observations of the bright radio sources J2253+1608 and J0555+3948 were included to serve as fringe-finders and to correct large delay errors.

The observing bandwidth of 32 MHz per polarization was divided into 4 sub-bands of width 8 MHz. Both senses of polarization were recorded with 1-bit sampling. The data were correlated in Socorro, New Mexico, producing 16 channels of width 0.5 MHz from each sub-band, with an integration time of one second.

Calibration was performed with aips444The Astronomical Image Processing System (aips) is developed and distributed by the National Radio Astronomy Observatory (NRAO). using standard procedures summarized as follows. Obviously corrupted data, and data taken at elevations less than $`15\mathrm{°}`$, were discarded. Visibility amplitudes were calibrated using on-line measurements of antenna gains, system temperatures, and voltage offsets in the samplers. Large delay errors were removed by fringe-fitting the data from a 1.5-minute observation of J2253+1608 and applying the delay corrections to all the data. At 2.3 GHz, we corrected residual rates, delays, and phases by fringe-fitting the data from J0134–0931 directly; at the higher frequencies we fringe-fitted data from J0141–0928 and interpolated the solutions. After calibration, the data were averaged in time and frequency in order to reduce the data volume, but we ensured that the final sampling was still fine enough to produce $`<`$1% smearing over the desired field of view.

### 3.2 Images and general results

Imaging was performed with aips. For each frequency, we first applied a Gaussian tapering function to the visibility data, to emphasize the low-resolution (short-baseline) information, and used the clean algorithm to create a preliminary model of the detected radio structures. This model was used to self-calibrate the antenna phases. A new map was created with no tapering function, and the process of cleaning and self-calibration was repeated until no further improvement was noted. We then used the latest clean model to self-calibrate the relative amplitudes of the antennas with a solution interval of 30 minutes. The amplitude corrections were typically smaller than 10%. We used a moderately uniform weighting (robust $`=0`$, in aips) for the final images.

We also applied the self-calibration and imaging procedures described above to the previously published 1.7 GHz and 5.0 GHz data, for the sake of uniformity. This resulted in an improvement in dynamic range over the previously published images, due to the increased care in removing corrupted data, the use of amplitude self-calibration, and the different choice of visibility weighting.

In all cases except 15 GHz, we created one wide-field image, centered between components B and D in right ascension, and between D and E in declination. At 15 GHz we used the multiple-field clean algorithm to deconvolve three fields simultaneously: one field centered on components A, B, and C; one field centered on component E; and one field centered on components D and F.

Figure 3 shows wide-field grayscale images for all frequencies except 15 GHz. These images provide an overview of the detected radio structures. To ease the visual interpretation of these images, they were created after applying a Gaussian tapering function to the visibility data, and were restored with a circular Gaussian beam. For the 1.7 GHz and 2.3 GHz images, the tapering function declined to 30% at a visibility radius of $`10^6\lambda `$, and the restoring beam had a diameter of 20 milli-arc seconds (mas). For the 5.0 GHz and 8.4 GHz images, the corresponding values were $`3\times 10^6\lambda `$ and 7 mas.

We also present narrow-field contour maps of the fields surrounding components A, B, and C (Figure 4), components D and F (Figure 5), and component E (Figure 6). The 15 GHz contour maps are presented separately for each component (Figure 7). All the contour maps were created without any tapering function and therefore employ the full resolution of the observations. These contour maps allow for a closer examination of individual component sizes and shapes. The most important parameters describing these full-resolution images are given in Table 1: the image dimensions, the major and minor axes (full-width-half-maximum) and orientation of the restoring beam, and the root-mean-squared (R.M.S.) noise level.

We draw attention to some features of these images and contour maps that will be discussed in detail, and placed in the context of the KW03 models, in subsequent sections:

1. Component F is detected at all frequencies except 15 GHz. It is especially prominent at 1.7 GHz but fades relative to D as the frequency increases (§ 3.3).
2. Components C and E are heavily resolved at 1.7 GHz, are faint and resolved at 2.3 GHz and 5 GHz, and are compact and bright at the highest frequencies. By contrast, components B and D are compact at all frequencies (§ 3.4).
3. In the 1.7 GHz map, there is extended radio emission near the expected location of C, but its centroid appears to be shifted $``$30 mas south (§ 3.5).
4. The radio arc between A and B was detected at the two lowest frequencies. In addition, at 1.7 GHz, there is a significant extension of component A towards the expected location of C (§ 3.6).
5. At 15 GHz (Fig. 7), component B appears to be a fainter and parity-reversed copy of component A. Component C is elongated along a line pointing to A (§ 3.6).

### 3.3 Component F

We measured the total flux density of each radio component by fitting an elliptical Gaussian function to the appropriate region of each VLBA image using the aips procedure jmfit, which takes into account the ellipticity of the restoring beam. The results are given in Table 2. The flux densities of components A, B, D, and F are plotted in Figure 8 as a function of frequency. Also plotted is the total flux density of all components, from the measurements compiled by W02. The upper limit on the flux of F at 15 GHz represents the $`3\sigma `$ level in that region of the image. Results for components C and E are not plotted because, for these components, a substantial amount of flux is “resolved out” at low frequencies due to the interferometer’s lack of short baselines.

The continuum spectra of components A, B, and D are nearly identical, while the spectrum for component F is significantly steeper. To make the difference in spectral slope obvious to the eye, we have divided the flux densities of D by 9, to match the 1.7 GHz flux density of F, and plotted the results with open squares. Using the measurements at which the components are fairly compact and the fitting results are therefore most accurate, we find spectral slopes of $`\alpha _\mathrm{A}=0.64\pm 0.08`$, $`\alpha _\mathrm{B}=0.65\pm 0.07`$, and $`\alpha _\mathrm{D}=0.63\pm 0.07`$ (where $`S_\nu \nu ^\alpha `$). These values agree with the value $`\alpha =0.69\pm 0.04`$ measured by W02 using lower-resolution data. For component F, using measurements at all frequencies for which it was detected, $`\alpha _\mathrm{F}=1.07\pm 0.13`$.

In the KW03 scenario, A–E are multiple images of the same quasar, and should have nearly identical spectral indices, but F need not have (and probably would not have) the same spectral index. It was already established by W02 that A–E have the same spectral index. By establishing that F has a different spectral index than the other components, our analysis has confirmed this key expectation of the KW03 models.

### 3.4 Scatter broadening

One objection to the KW03 models, or any scenario in which components A–E are gravitationally lensed images of a single source component, is that components C and E were observed to have a much lower surface brightness at 1.7 GHz and 5.0 GHz than the other components. The problem is that gravitational lensing cannot alter surface brightness. However, it is possible for a non-gravitational effect, namely scatter broadening due to electron-density fluctuations in the lens galaxy or our own Galaxy, to cause differences in the surface brightness of lensed images. If different images pass through different electron columns they will be broadened by different amounts.

Scatter broadening can be recognized because the angular size of an affected source increases strongly with observing wavelength. For a source that is intrinsically a point source, the angular size due to scatter broadening is typically proportional to $`\lambda ^n`$ where $`n2`$ (for a review, see Rickett, 1990). A particularly well-studied case is Sgr A (see, e.g., Jauncey et al., 1989; Lo et al., 1993; Doeleman et al., 2001). Among gravitational lenses, the best-studied case of differential scatter broadening is PKS 1830–211 (Jones et al., 1996; Guirado et al., 1999). In that case, the broader image passes through the spiral arm of the lens galaxy where a higher electron column density might be expected (Winn et al., 2002b; Courbin et al., 2002).

Our VLBA maps strongly suggest that components A, C, and E are being scatter broadened. All three of these components are fairly compact at 15 GHz, but are heavily resolved at 5 GHz and very heavily resolved at 1.7 GHz. This is despite the decreasing angular resolution of the observations as one proceeds to lower frequencies. To test the $`\lambda ^2`$-dependence, we fitted elliptical Gaussians to all the components detected in each image (as in § 3.3). The fitting procedure takes into account the intrinsic ellipticity of the beam and finds the deconvolved major and minor axes ($`b_{\mathrm{maj}}`$ and $`b_{\mathrm{min}}`$). We note that it would be preferable to fit models to the visibility data rather than the image data, because the noise properties of the visibility data are better understood, but in this case visibility-fitting was unfeasible because of the large data volume. We also note that not all the components are well described by ellipses, especially at 15 GHz, but the results from elliptical fits should be accurate enough to determine the approximate logarithmic scaling.

Figure 9 shows the dependence of $`b_{\mathrm{maj}}`$ on $`\lambda `$ for components A, C, and E. As a measure of the angular resolution, we have plotted black dots indicating the length of the minor axis of the restoring beam, which is approximately proportional to $`\lambda `$. For comparison, we have also plotted the results for component D, which does not exhibit scatter broadening; its size is comparable to or smaller than the beam minor axis in all cases except 1.7 GHz. (Component B does not appear to be scatter-broadened either, but the results are not plotted in order to reduce clutter.)

Components A, C, and E are well-approximated by the $`\lambda ^2`$ law (dashed line), although the true dependence is not as steep, especially when the 15 GHz (2 cm) point is included. Given that scatter broadening seems inevitable as an explanation for the wavelength-dependent morphology of these components, we believe that the departure from the $`\lambda ^2`$ law shows that the observed angular size is partially intrinsic at the shortest wavelengths. The dotted line shows the expectation of a scatter-broadened source with intrinsic size 2 mas. The source shapes in the 8.4 GHz and 15 GHz images appear to be mainly intrinsic. The flux densities, major and minor axis sizes, and position angles of the components, taken from the 15 GHz image in order to minimize scatter broadening, are given in Table 3.

The observed angular sizes can be used to estimate the scattering measure (SM) along the lines of sight to the broadened components. The SM is defined as the line-of-sight integral of $`C_n^2`$, where $`C_n^2`$ is the normalization of the power spectrum of electron-density fluctuations. Assuming a Kolmogorov spectrum, the SM can be related to the observed broadening of an extragalactic source via the relation (see Fey, Spangler, & Cordes 1991, Eq. 4; or Cordes & Lazio 2003, Eq. A17):

$$\mathrm{SM}=\left(\frac{\theta _{\mathrm{FWHM}}}{130\mathrm{mas}}\right)^{5/3}\left(\frac{\nu }{1\mathrm{GHz}}\right)^{11/3},$$
(1)
where $`\nu `$ is the frequency at the lens redshift, and SM has the units m-20/3 kpc. This predicts that angular size should vary as $`\lambda ^{2.2}`$, a little steeper than observed. Nevertheless, the best fit to the data plotted in Fig. 9 gives SM $``$ 15–25 m-20/3 kpc. This value is too large to be attributed to ionized material in our Galaxy. The Galactic latitude of J0134–0931 is $`70\mathrm{°}`$, and the SM along such high-latitude lines of sight in the Galaxy is typically less than $`10^3`$ m-20/3 kpc (Cordes & Lazio 2003, Fig. 6).

By contrast, the SM for low-latitude lines of sight in our Galaxy can exceed $`10^3`$ m-20/3 kpc, and has an average value of order unity within $`|b|<5\mathrm{°}`$. This suggests that a spiral lens galaxy viewed nearly edge-on could produce the observed broadening, although it seems unlikely that the disk would be aligned just right to intersect all three components A, C, and E. It could also be that the lens galaxy has a larger electron content than our Galaxy. We caution, however, that comparisons with our Galaxy may be misleading because the observed major axes of $``$30 mas at 18 cm correspond to a physical size of $``$150$`h^1`$ pc at the lens redshift, which is much larger than the length scales probed by scattering effects in our Galaxy. (For exactly that reason, observations of scatter broadening in gravitational lens systems such as J0134–0931 may prove useful for characterizing electron-density fluctuations on scales inaccessible in our Galaxy.)

### 3.5 Counter-images of component F

Component F, in the KW03 models, is actually the brightest of three images of the second source component S2. The other two images of S2 (the “counter-images” of F) are expected to be located just south of components C and E. For convenience, we refer to the hypothetical counter-images near C and E as components $`\mathrm{C}^{}`$ and $`\mathrm{E}^{}`$, respectively. The KW03 models suggest that the flux ordering, from brightest to faintest, is F–$`\mathrm{C}^{}`$$`\mathrm{E}^{}`$. It is also possible that there are two additional counter-images located between A and B, for a total of five, but the models that produce these two additional counter-images occupy a tiny volume in the phase space of model parameters explored by KW03.

There is suggestive evidence for one of these counter-images in the 1.7 GHz map. In Fig. 4, black dots mark the positions of A, B, and C (as measured at higher frequencies), and the predicted position of $`\mathrm{C}^{}`$. There is a puff of extended radio emission located $``$30 mas due south of the expected position of component C. The location of the peak emission is not consistent with C, but it is within 6 mas of the predicted position of $`\mathrm{C}^{}`$. Its total flux density is $``$10 mJy, or $``$90% of the total flux density of F, which is also the expectation of the KW03 models.

It is therefore possible that C is almost completely resolved out at 1.7 GHz, and the extended emission that is observed is instead the scatter-broadened counter-image $`\mathrm{C}^{}`$. The reason we consider the evidence “suggestive” rather than definitive is that $`\mathrm{C}^{}`$ was not detected at any other frequency. In particular, the 8.4 GHz data are inconsistent with $`S_\mathrm{C}^{}/S_\mathrm{F}=0.9`$, unless $`\mathrm{C}^{}`$ is being scatter-broadened even at that higher frequency. Therefore, we withhold final judgment on whether the detected flux is actually $`\mathrm{C}^{}`$ until it can be detected at a higher frequency. The challenge will be to achieve the required signal-to-noise ratio in spite of the steep spectrum of S2 (see § 3.3).

Assuming that we did not detect additional radio components, we can set upper limits on the fluxes on $`\mathrm{C}^{}`$ and $`\mathrm{E}^{}`$. We express these limits as ratios relative to component F, because lens models directly prescribe the magnification ratios of the counter-images relative to component F. We derive the limits from the 8.4 GHz map, because that is the highest-frequency image in which F is detected. By requiring the flux density of a point source be below $`5\sigma `$, where $`\sigma `$ is the noise level of the 8.4 GHz image, we find $`S_\mathrm{C}^{}/S_\mathrm{F}<0.57`$ and $`S_\mathrm{E}^{}/S_\mathrm{F}<0.57`$.

As discussed in detail by KW03, the limit on $`\mathrm{E}^{}`$ is not constraining, but the limit on $`\mathrm{C}^{}`$ formally rules out all but a few models. However, as mentioned previously, the true limit on $`\mathrm{C}^{}`$ may be weaker due to scatter broadening. For example, assuming that the observed sizes of C and E at 8.4 GHz are due entirely to scatter broadening, we might allow $`\mathrm{C}^{}`$ and $`\mathrm{E}^{}`$ to have the same ratio of total to peak flux density. Modified in this way, the upper limits become $`S_\mathrm{C}^{}/S_\mathrm{F}<2.1`$ and $`S_\mathrm{E}^{}/S_\mathrm{F}<1.1`$, which are consistent with all the KW03 models. The truth lies somewhere between these extremes, because the observed sizes of C and E at 8.4 GHz are probably at least partially intrinsic (see § 3.4).

### 3.6 Extended emission

The radio arc joining components A and B was known previously, and is a natural consequence of the KW03 models. In addition, the KW03 models predict there should be a bridge of radio emission between components A and C. The detailed surface brightness profile of the bridge depends on the source size and shape, but it is otherwise a robust prediction.

There are two indications in the radio maps that this is correct. The first indication is found in the lowest-frequency (1.7 GHz) image, in which there are faint traces of extended flux between components A and C. This can be seen in the upper left panel of Fig. 3, as a faint prong extending eastward from component A, but is more easily seen in the upper left panel of Fig. 4. The lowest 3 contours ($`3\sigma `$, $`6\sigma `$, and $`12\sigma `$) bulge eastward from A to C.

The second indication is the orientation of component C in the higher-frequency (5 GHz, 8.4 GHz and 15 GHz) images. The deconvolved position angle of C is oriented nearly along the direction to component A. For example, at 15 GHz, where we expect the observed orientation to be intrinsic, component C is elongated with position angle $`100\mathrm{°}\pm 6\mathrm{°}`$, whereas the direction from C to A has position angle $`110\mathrm{°}`$ (see the lower left panel of Figure 4).

### 3.7 Component shapes

Since we have argued that the observed sizes of the radio components at 15 GHz are mainly intrinsic, the morphologies of the components should provide clues about the properties of the background source and the correct lens model. In Fig. 7, component A appears to be a two-sided core-jet source. There is a longer and more prominent jet directed southwest of the core, towards B, and there is a shorter jet directed northeast. Component B also appears to be a core-jet source, but only one jet is observed, and is directed northeast, towards A. Apparently, A and B are opposite-parity images of a two-sided core-jet source, but only one jet is apparent in B due to the relative compression of that image by the lens mapping.

The morphologies of the other components are simpler. Component C is highly elongated towards A, as already noted. Component D is mainly compact but does have a faint extension to the southwest, suggesting that D and A have the same parity. Component E is resolved but is too faint for a detailed examination. Component F was not detected.

The KW03 models predict that A and B have opposite parities, and that A and D have the same parity, which is consistent with the observed morphologies. In addition, the position angles given in Table 3 for components A–E can be produced by 23% of the models discussed by KW03. However, we caution against over-interpreting the shape parameters. As described in § 3.4, they were derived by the simple procedure of fitting a single elliptical Gaussian function to each component. The quoted errors are the statistical errors in the fitting process, and are therefore under-estimated for those components that are not well described by ellipses (especially A and B). The most secure parameters are the position angles, and the least secure parameters are the minor axis lengths (which in some cases are considerably smaller than the restoring beam).

## 4 Optical evidence

### 4.1 Observations and data reduction

On 2001 July 18 we obtained optical images of J0134–0931 in the F555W (hereafter, $`V`$) and F814W ($`I`$) filters on the WFPC2 camera of the HST. We obtained two 21.7-minute exposures at each of two dither positions, for a total exposure time of 86.8 minutes per filter. The target was centered in the PC chip, which has a pixel scale of $`0\stackrel{}{\mathrm{.}}0456`$ pixel-1. The exposures were combined and cosmic rays rejected using standard iraf555The Image Reduction and Analysis Facility (iraf) is a software package developed and distributed by the National Optical Astronomical Observatory, which is operated by aura, under cooperative agreement with the National Science Foundation. procedures. The images are shown in the upper two panels of Figure 10. Components A and B were detected in both filters. Component D was detected in $`I`$-band but is nearly absent in $`V`$-band.

These same data were previously analyzed by Hall et al. (2002), who also measured the optical flux ratios and detected residual flux after subtracting point sources at the positions of the radio components. Our results confirm this previous work. However, as described below, our more detailed methods have allowed us to proceed further in the interpretation of the results.

### 4.2 Detection of the lens galaxies

We built photometric models of the HST images using custom software for convolving a parametric model with a simulated PSF and varying the model parameters to achieve the best fit to the images. The simulated PSF was generated with the Tiny Tim program by Krist & Hook (1997). For the $`V`$-band image, the model consisted of 5 point sources, representing the radio components A–E. At $`I`$-band we found that components A and B are not well represented by point sources, possibly due to a contribution from the quasar host galaxy. We obtained much smaller residuals by allowing for small Gaussian components at the locations of both A and B, in addition to point sources.

Even with this modification, there were still significant residuals. In Figure 10 we show the residuals as grayscale images overlaid with white contours. The positions of the subtracted components are marked with white circles. The ring of negative residuals surrounding component A is spurious and due to imperfect PSF subtraction, but there are two peaks of positive emission in the residual images.

1. The “northern peak,” $`0\stackrel{}{\mathrm{.}}14`$ south and $`0\stackrel{}{\mathrm{.}}12`$ east of E, and
2. The “southern peak,” $`0\stackrel{}{\mathrm{.}}37`$ south and $`0\stackrel{}{\mathrm{.}}03`$ east of C.

We believe these two peaks represent genuine emission rather than artifacts of the subtraction. This pattern of residuals did not appear when we applied the same fitting and subtraction procedure to an isolated field star on the PC chip. In addition, although the peaks are more prominent in $`I`$-band, they are detected in both filters, despite the completely different PSF.

Hall et al. (2002) already noted the presence of low-surface-brightness flux underlying the point sources. However, the pattern in their residual image is less obvious, possibly due to an inferior PSF subtraction. They used only integral pixel shifts, did not simultaneously fit for the fluxes of the point sources, and did not allow for any extension of A and B.

Interestingly, G02 also found two residual peaks of emission by deconvolving a ground-based $`K^{}`$-band image. The locations of each $`K^{}`$-band peak is within $`0\stackrel{}{\mathrm{.}}2`$ of the corresponding peak in our images, which is fairly good agreement, considering that the spurious residuals around A may be corrupting our measurements of the peak locations. Although G02 were unsure whether the $`K^{}`$-band peaks represented genuine emission, or artifacts of the deconvolution, we argue that the observation of similar peaks at both $`I`$-band and $`V`$-band lends confidence to all three detections.

The detection of two diffuse light sources is a further confirmation of the two-galaxy KW03 lens models. As discussed further by KW03, the predicted positions of the galaxies are both within $`0\stackrel{}{\mathrm{.}}15`$ of the observed locations at $`I`$-band and $`V`$-band. The approximate range of lens galaxy positions predicted by KW03 is indicated by the two small gray ellipses in Figure 10.

It is not clear how seriously to take the observed morphologies of the northern and southern peaks, because they are faint, and because of possible contamination by spurious residuals from the PSF subtraction. Nevertheless, given the morphology of the southern peak observed in both filters, we conjecture that the southern lens galaxy is elongated east-west, and overlaps component D at its easternmost extent. This would explain the very red color of D (see § 4.3).

### 4.3 Photometry and colors of the components

We measured the optical fluxes of components A–E using the PSF-fitting procedure described above. The optical flux ratios for A–E, relative to component B, are given in Table 4 for both the $`I`$ and $`V`$ filters. For some components we could derive only upper limits. For reference, we also list the radio values measured by W02.

To convert the fluxes to standard photometric magnitudes, we adopted zero points of 21.69 for F814W and 22.56 for F555W, based on the work of Dolphin (2000) and corrected for infinite aperture. No correction was attempted for the charge-transfer-efficiency (CTE) problem of WFPC2. The total magnitudes of all the components in our model was $`I=18.84\pm 0.05`$ and $`V=22.54\pm 0.07`$, in good agreement with the ground-based photometry of W02 ($`I=18.79\pm 0.05`$, $`V=22.47\pm 0.09`$).

All the components that were detected are very red. Even component B, the bluest of the detected components, has $`VI=3.4`$. This is redder, for example, than any of the 157 quasars in the radio-selected (and optically unbiased) sample of Parkes sources that was investigated by Francis, Whiting, & Webster (2000) (see Fig. 4 of that work).

The optical flux ratios disagree strongly with the radio flux ratios, even though gravitational lensing is achromatic. This is a common observation among gravitational lenses, usually attributed to dust extinction by the lens galaxy (see, e.g., Falco et al., 1999) or optical microlensing (Wambsganss & Paczynski, 1991), although in this case the differential extinction is larger than usual. By comparing the radio flux ratios and the optical flux ratios, we can estimate the difference between the $`I`$-band extinction $`A_I`$ incurred by B (the bluest component), and the $`I`$-band extinction incurred by the other components:

$$\mathrm{\Delta }A_I(\mathrm{A})=A_I(\mathrm{A})A_I(\mathrm{B})=2.5\mathrm{log}_{10}\left(\frac{f_I(\mathrm{A})/f_I(\mathrm{B})}{f_{\mathrm{radio}}(\mathrm{A})/f_{\mathrm{radio}}(\mathrm{B})}\right),$$
(2)
where $`f_I(\mathrm{B})`$ refers to the $`I`$-band flux of component B, and so forth. The results for the components, in order from bluest to reddest, are: $`\mathrm{\Delta }A_I(\mathrm{B})0`$, $`\mathrm{\Delta }A_I(\mathrm{A})=0.97`$, $`\mathrm{\Delta }A_I(\mathrm{C})>1.7`$, $`\mathrm{\Delta }A_I(\mathrm{D})=3.1`$, and $`\mathrm{\Delta }A_I(\mathrm{E})>4.4`$. (The position of C in the ordering is uncertain because the result is a relatively weak upper limit.)

The most plausible explanation for the different optical colors of all the lensed components is differential reddening due to dust in the lens galaxies. The main alternative explanation, microlensing, is implausible because it is unlikely that microlensing would affect more than one image to the required degree at the same time. Wyithe & Turner (2002), for example, estimate there is only a 30% chance for even one image of a multiple-image quasar to vary by more than 0.5 mag over 10 years. In addition, the dust hypothesis explains why the fainter components are also redder; this would be a coincidence under the microlensing hypothesis. Therefore, at least part of the reason for the extreme red color of J0134–0931 is dust extinction in the lens galaxies. This need not be the entire reason. Because we can determine only the differential extinction, rather than the total extinction, we cannot rule out the existence of dust in the host galaxy that reddens the background source, or an intrinsically red background source.

## 5 Summary and conclusions

The evidence presented in this paper strongly supports the models developed in detail by Keeton & Winn (2003): the unique radio morphology of J0134–0931 represents five images (A–E) of a single radio-loud quasar. There is also a sixth radio component (F) representing a different part of the background radio source.

Radio observations with the VLBA at five widely spaced frequencies demonstrate that F has a steeper spectral index than components A–E, ruling out six-image models and supporting the notion that F represents a different part of the background radio source. One of the predicted counter-images of F may have been detected at 1.7 GHz, but this needs to be confirmed with more sensitive observations at higher frequencies. If there remains any doubt about the KW03 scenario, detection of both of the counter-images would be definitive.

The multi-frequency maps also demonstrate that the lower surface brightness of components A, C, and E at long wavelengths is due to differential scatter broadening. This defends the theory that they are multiple images of a single quasar from the objection that surface brightness must be conserved by gravitational lensing. Components B, D, and F do not appear to be significantly scatter broadened.

At short wavelengths, where scatter broadening is minimized, the resolved shapes of the radio components can be broadly explained by a subset of the KW03 models. In particular, component C is aligned in the direction of component A, as predicted generically by those models. The radio arc between components A and B, which was previously detected at 1.7 GHz and which we also detected at 2.3 GHz, emerges naturally from the models. At the longest wavelength there is a trace of extended emission between A and C, again fulfilling the expectation of the lens models that these components should be connected by an arc.

The positions of the two lens galaxies proposed by KW03 are in rough agreement with two residual peaks in optical images with the HST, and also with two residual peaks identified by G02 in a deconvolved $`K^{}`$-band image. The northern galaxy appears to be responsible for the reddening and scatter broadening observed in components A, C, and E. The southern galaxy appears to be responsible for the reddening of component D. Component B, which is the most distant in projection from either galaxy, is not scatter broadened and is the bluest component, as expected. It should be possible to improve the characterization of the lens galaxies with space-based images, or ground-based images with a very high signal-to-noise ratio (allowing a reliable deconvolution). Given the presence of dust and ionized material in the lens galaxies and the possibly elongated morphology observed for the southern lens galaxy, it seems likely that both galaxies are spiral galaxies.

We thank Ramesh Narayan and David Rusin for useful discussions. We are indebted to Steven Beckwith for providing Director’s Discretionary Time on the HST for this project, and we thank Emilio Falco, Chris Impey, Brian McLeod, and Hans-Walter Rix, for helping to obtain the HST data as members of the castles collaboration. J.N.W. is supported by an NSF Astronomy & Astrophysics Postdoctoral Fellowship, under grant AST-0104347. C.S.K. is supported by nasa ATP grant NAG5-9265 and HST grant GO-9133.


format: text

text:
# Core dissolution and the dynamics of massive stars in young stellar clusters

## 1 Introduction

The origin of massive stars is a poorly understood phenomena. One clue to their formation is that they are generally found in the dense cores of young stellar clusters (Hillenbrand 1997; Carpenter & Sanders 1998; Clarke et al. 2000; Mermilliod & García 2001; Zinnecker et al. 1993). Their location in these very young systems cannot generally be explained by dynamical mass segregation (Bonnell & Davies 1998). Higher than observed primordial core densities (Kroupa 2002; Kroupa et al. 2001) can cause rapid mass segregation in the core region. The central location of massive therefore must tell us something about their birth. One possibility is that the accretion rates in the centre of the cluster are higher due to the deeper potential wells (Bonnell et al. 2001; Bonnell & Bate 2002) and that therefore the protostars that are present in the core are able to accrete up to a high mass (eg. Behrend & Maeder 2001, Klessen 2001). An additional possibility is that as the stars in the core accrete, the core contracts due to the increase in the potential energy. The increase in stellar density may then be sufficient such that collisions between intermediate mass stars occur to form the massive stars (Bonnell, Bate, & Zinnecker 1998; Bonnell & Bate 2002). One of the successes of these models is that they naturally account for the formation of massive stars in the centres of stellar clusters. These models predict the location of massive stars at the end of the accretion, mass-buildup phase. At this point the cluster contains significant mass in gas which can still affect the cluster dynamics.

Once a massive star has formed in a cluster, it can have dramatic effects on its surroundings. The strong winds and ionising photons inject significant energy into the surrounding gas which will begin to expand, eventually being cleared from the cluster. This gas removal has two effects: firstly, it halts any residual accretion or star formation; secondly, if the gas still comprises a significant fraction of the total mass, the cluster will expand to reflect the new, lower potential. Simulations of this process (Geyer & Burkert 2001; Kroupa et al. 2001; Boily & Kroupa 2003) have concentrated on the global evolution of the cluster as the gas is removed. They find that if the gas is a small fraction of the total mass, or if it is removed slowly relative to the dynamical timescale, then the cluster can retain significant fraction of its initial stars. Kroupa et al. (2001) suggests that early dynamical mass segregation of the most massive stars can survive during and after a gas expulsion phase. Kroupa (2000) has presented results specifically for models of the Orion Nebula Cluster (ONC) including the gas expulsion scenario reflecting the situation we model here, but without primordial mass segregation.

In this paper, we are concerned with the dynamical evolution of the cluster core when the gas is removed. Gas removal occurs from inside out as the massive stars inject energy into their environment. The core of the cluster will be cleared first and the core stars will be affected by the gas removal before the rest of the cluster. This has serious implications for the dynamics of the massive stars, their location in the cluster, or even their possible ejection from the cluster. We address these questions using numerical simulations presented below.

## 2 Gas expulsion from core

Gas removal occurs due to the injection of energy from the massive stars in the core. The energy produced by these O and B type stars can take the form of strong stellar winds of up to a few $`\times 1000`$ km s-1 (Churchwell 1999) or the formation of an expanding HII region ($`v_{HII}10`$ km s-1). In both cases, the expansion velocity is significantly higher than the velocity dispersion seen in young stellar clusters due to the higher initial density. If the core of the cluster is about 0.1pc this means that the crossing time for the core $`(t_{\mathrm{cross}}5\times 10^4yr)`$ is at least an order of magnitude larger than the gas expulsion time. We may then assume that the period of gas expulsion has negligible effect on the dynamics of the constituent stars in the core, and so we assume instantaneous gas expulsion. Thus we incorporate an initial surplus of kinetic energy into the core stars which results in their being no longer in virial equilibrium. This implies that the gas is only removed from the core as stars further out remain in virial equilibrium. The timescales $`(10^5yr)`$ and radial extent $`(0.1pc)`$ reflect the situation in an ultra-compact HII region where a stellar wind has evacuated a small inner region before stalling due to mass loading (Lizano et al. 1996), ram pressure of the infalling gas (Wood & Churchwell 1989), or other possible mechanisms (Churchwell 1999).

## 3 Numerical Simulations

Over the timescales involved it is only gravitational forces which play a significant rôle in the subsequent evolution. Thus the system is sufficiently modelled with a pure N-body simulation. The numerical simulations reported here were performed using the NBODY2 code (Aarseth 2001). This code is extremely efficient and accurate at following a stellar dynamical system for tens to a hundred dynamical times. It includes the AC neighbour scheme (Ahmad & Cohen 1973) and a softened potential (here $`r_{\mathrm{soft}}0.1\times r_{\mathrm{core}}`$) which minimises computational expense at the cost of not following any (close) binary systems. As we are following the stellar dynamics over only a few local dynamical times where distant interactions dominate, NBODY2 is sufficient for our needs.

Our cluster is assumed to be in a virialised state prior to gas expulsion, from arguments of the star formation timescale compared to the crossing times, typical systems may indeed be near to viral equilibrium (Boily & Kroupa 2003). The stellar system consists of 1500 stars populated with masses distributed according to a Salpeter mass function, $`\varphi (m)m^{2.3}`$, ranging from $`0.110\mathrm{M}_{}`$ and distributed in phase space with a plummer density profile. Additionally, it is necessary to model an initially mass segregated system such that the massive stars are found in the core region of the cluster, as described above. We construct such initial conditions by allowing our randomly distributed stellar system to relax for several dynamical times until the massive stars have sunk to the core. This system is then used as our initial condition for the simulations in this investigation. It must be emphasised that this process of achieving initial conditions is only to construct a model of the observed initial state and not to represent the physical method by which massive stars end up in the cores of very young stellar clusters in the first place. Our pre-initial clusters are evolved until little more mass segregation can occur. The resulting initial mass distribution in the core is shown in Figure 1 showing that the high mass stars have sunk to the core.

Our results are presented scaled to the physical units which apply to the Orion nebula cluster as taken from Hillenbrand & Hartmann (1998), it should be emphasised that our models are not intended to model the ONC which has a much wider mass range $`(0.150\mathrm{M}_{})`$. The half mass radius of the cluster, $`r_{\mathrm{half}}=0.8\mathrm{pc}`$. We choose the core radius to be $`r_{\mathrm{core}}=0.08\mathrm{pc}`$ being 0.1 the half mass radius, it is this region within $`r_{\mathrm{core}}`$ in which we consider the gas expulsion to be effective. The mass of our cluster is $`1200\mathrm{M}_{}`$, and N$`=1500`$ so that the mean stellar mass, $`<\mathrm{M}_{}>=0.8\mathrm{M}_{}`$. Finally, we follow the evolution of the core stars to $`t1.0\mathrm{Myr}`$ $`(20t_{\mathrm{cross}})`$, by which time the dynamical evolution we are interested in will have occurred.

## 4 Run Parameters

The effect of gas removal is included by way of a super-virial distribution of the stellar velocities inside the cluster’s core. The kinetic energy of these stars is taken to be that required to balance an additional mass of gas, $`M_{gas}`$. We have parametrised our model by the ratio of $`M_{gas}/M_{stars}`$ prior to gas expulsion in the core. The values of this parameter used and the equivalent star formation efficiencies are shown in Table 1. This method mimics the case when gas removal is instantaneous. This is the worst case scenario for retaining the stars and thus gives a conservative estimate of how many massive stars will be retained in the core.

In practice we take our initial system as described above and increase the kinetic energies of the stars in the core by the following energy argument,

$$[M_{gas}+M_{stars}]v_{old}^2=M_{stars}v_{new}^2$$
(1)
we can then directly scale the core velocities such that,

$$v_{new}^2=\frac{v_{old}^2}{ϵ}$$
(2)
where we have introduced the star formation efficiency, in terms of the gas ratio as,

$$ϵ=1/(1+\frac{M_{gas}}{M_{stars}}).$$
(3)
For each value of the initial parameter, $`M_{gas}/M_{stars}`$, 8 initial realisations were constructed using a different random seed for each. This helps to extract the mean behaviour of the massive stars in the clusters and not individual fluctuations due to their small numbers.

## 5 Results

The purpose of our investigation is to discover whether a mass segregated core will survive after a gas expulsion phase, and also to gain general insight to the subsequent motions of the massive stars. The following results examine properties of the core: mass segregation, average mass within the core, fraction of stars remaining in the core, and core survivability. We also look at the phase space evolution of the individual massive stars by examining their radial evolution and their angular momentum.

### 5.1 Mass Segregation

We investigate the evolution of mass segregation in the core, which, in this case, amounts to the survival of the most massive stars in the core.

After $`0.62\mathrm{Myr}`$ $`(12t_{\mathrm{cross}})`$ we take the ratio of the number of stars within $`r`$ to the number initially contained within given radii ($`R_{\mathrm{core}}`$ and $`2R_{\mathrm{core}}`$). This cumulative number ratio is calculated independently for 3 stellar mass ranges: $`M_{}1.08\mathrm{M}_{},\mathrm{\hspace{0.33em}1.08}\mathrm{M}_{}<M_{}<5.4\mathrm{M}_{},\mathrm{\hspace{0.33em}5.4}\mathrm{M}_{}M_{}`$, shown as a thin line, a medium line, and a thick line respectively in Figures 2(a)–(d). These figures show the cumulative number ratio for a $`M_{gas}/M_{stars}`$ parameter of (a) 0.25, (b) 1.0, (c) 3.0, and (d) 5.0. The top and bottom panels in each figure show the cumulative number ratio normalised to $`2R_{\mathrm{core}}`$ and $`R_{\mathrm{core}}`$ respectively. The initial state of the system at $`t=0`$ is shown as a dashed line in each figure. These figures are the combined results of all realisations for each particular $`M_{gas}/M_{stars}`$.

In particular, we are interested in the highest mass range, which we consider to represent the OB type stars in the cluster. Clusters with lower values of $`M_{gas}/M_{stars}`$ (figures 2(a) and (b)) show little significant change in the distribution of the massive stars, in fact the lower mass stars decrease in number more than the massive stars, making the core slightly more segregated. However, as soon as the effective kick to the core stars becomes large enough to boot out the massive stars we see a dramatic decrease in all the stellar numbers. Figures 2(c) and (d) show a decrease in the actual massive star numbers by 20% and 70% respectively.

An alternative and perhaps more intuitive way of looking at the mass segregation is to look at the average stellar mass within a certain radius, or as we do in this case the 100 stars nearest to the centre of mass. Figure 3 shows there is an almost linear decrease in $`\overline{M_{100}}`$ for each value of the $`M_{gas}/M_{stars}`$ parameter after $`M_{gas}/M_{stars}=1.0`$. Also shown, as a vertical bar on each point, is the range due to the variation in the initial realisations for each parameter. This shows clearly that as the stars are ejected out into the rest of the cluster, the massive stars are less of a characteristic of the central cluster regions. Although the mean stellar mass is slightly greater than the cluster mean $`<M_{}>=0.8\mathrm{M}_{}`$, at $`M_{gas}/M_{stars}=5.0`$ the mass segregation is much less significant.

### 5.2 Radial evolution of Massive stars

The question now arises of what becomes of the massive stars which the core has ejected. For comparison with observed clusters we would like an idea of where the stars are located in our simulations. We examine the radial distances of the most massive stars from the centre of mass with respect to time. Typical examples are shown in Figure 4.

These figures illustrate typical evolution of the 10 most massive stars which were initially in the cluster core, each figure showing one realisation for each value of the parameter $`M_{gas}/M_{stars}`$. As we have seen in the previous section, for a larger $`M_{gas}/M_{stars}`$, a greater number of massive stars leave the core. Figure 4 shows this trend explicitly and we can see the subsequent radial motion of the stars after gas expulsion. The grey region in the figures is the initial core in which the gas is assumed to have been removed. The half mass radius, $`R_{\mathrm{half}}`$, is also shown as a solid line at 0.8pc.

For a $`M_{gas}/M_{stars}=0.25`$, the figure shown top left is typical for all realisations and similar results are seen for $`M_{gas}/M_{stars}=1.0`$ (top middle). For $`M_{gas}/M_{stars}=2.0`$ (top right) it is more typical to see a couple of the stars ejected from the core but within the half mass radius. For $`M_{gas}/M_{stars}=3.0`$, 2 to 4 stars are ejected beyond $`R_{\mathrm{half}}`$ and another couple beyond $`R_{\mathrm{core}}`$ (bottom left). Figure 5 shows a resulting configuration in the $`xy`$ plane. There is a resemblance to the configuration of the massive stars in the Orion Nebula Cluster (ONC), showing that this is a possible mechanism leading to ONC type structures, where some of the massive stars are located a significant distance from the core. Although without observational velocities we cannot say how likely it is.

The tendency when $`M_{gas}/M_{stars}=4.0`$ (bottom middle) is for more massive stars to be completely ejected and still two or so to be found between $`R_{\mathrm{core}}`$ and $`R_{\mathrm{half}}`$. It is likely that most stars are ejected from the core for $`M_{gas}/M_{stars}=5.0`$ and also beyond $`R_{\mathrm{half}}`$ (Figure 4, bottom right).

### 5.3 Velocity Evolution

Observational measurements of velocities (proper motions as well as line-of-sight velocities) are now feasible with space borne instruments, such velocities may provide strong evidence to decide on the spatial origin of the stars. It is therefore of interest to explore the kinematics that result from gas expulsion and the ejection of massive stars from the core of the cluster.

The expectation is that stars which receive significant velocity increases escape the cluster on near-radial orbits. Instead, we show below that these stars generally have more tangential than radial motions. This can be explained by the interactions that occur due to the discrete nature of the stellar cluster. We can suppose that the orbits of these stars will encounter a great range of stars on highly varied orbits, its velocity signature will vary greatly, leaving the angular momentum measured with respect to the centre of the core density peak, $`|𝑳|`$, as the most robust general measure of the stars’ subsequent evolution. This measure of angular momentum is not affected by the motion of core over the duration of these simulations.

Figure 6 shows the ratio of the angular momentum at time $`t=0.62\mathrm{Myr}(12t_{\mathrm{cross}})`$ to the initial value immediately after the velocity kick $`|𝑳_T|/|𝑳_0|`$. All stars over $`5\mathrm{M}_{}`$, in all realisations of the parameter $`M_{gas}/M_{stars}`$ are shown. We have distinguished the stars by their radial position at time $`t=0.62\mathrm{Myr}`$: those within $`R_{\mathrm{core}}`$ are shown as small open circles ($``$), those beyond $`R_{\mathrm{core}}`$ shown as filled stars ($``$), and those beyond $`5R_{\mathrm{half}}`$ (those which have left the cluster) are shown as open triangles. Each radial range is shifted slightly to left or right for clarity.

The first effect to notice is the increasing spread of angular momentum with increasing initial $`M_{gas}/M_{stars}`$, such that for $`M_{gas}/M_{stars}=5.0`$ the gain in $`|𝑳|`$ for some stars reaches $`1000`$. The other striking feature is the bimodality between stars that have stayed in the core and stars which have left it. Stars remaining in the core have angular momenta distributed around their initial values, although plenty of angular momentum exchange has occurred resulting in $`100>|𝑳_T|/|𝑳_0|>0.01`$. The core escapers, however, almost exclusively increase their angular momenta such that $`1000>|𝑳_T|/|𝑳_0|>1`$.

Such a result indicates that the higher velocity massive stars must have encountered a greater number of stars with a broader range of orbits on their journey outward through the cluster. A combination of their gain in $`|𝑳|`$ and high velocity puts many stars on large tangential orbits which do not pass through the core. This was seen in Figure 4 as stars are ejected from the core but are seen to orbit beyond $`R_{\mathrm{core}}`$. As a consequence of this the velocities of the ejected massive stars will be seen to have predominantly tangential motions to the cluster core.

There is a further subset of stars which appear to leave the cluster entirely (a limit of $`5R_{\mathrm{half}}`$ is used to identify cluster leavers), these are shown as open triangles in Figure 6. The physical velocities of these stars range from $`10`$ to $`40\mathrm{km}\mathrm{s}^1`$. A result which allows us to rule out the existence of observed high velocity stars ($`200\mathrm{km}\mathrm{s}^1`$) as a consequence of the two body interactions which occur in this scenario.

### 5.4 Core Survival

To summarise the results in terms of the original question posed — how does gas expulsion from the core affect the core’s survival? — we have plotted (Fig 7) the fraction of massive stars remaining in the core for each value of $`M_{gas}/M_{stars}`$. It is perhaps a matter of definition as to what constitutes a surviving core, and in this case we take a value of 50% of the original massive stars. More precisely we say that if more than half the core (or rather its constituent massive stars) is dissipated subsequent to the gas expulsion phase then the core does not survive. This is supported by the fact that the actual number of massive stars in the core are of the order of a few, once half of these have left the system there may be only one or two left, not enough to constitute a core.

The relation between $`M_{gas}/M_{stars}`$ and the core survivability is fairly smooth (figure 7). The core is relatively undisturbed up to $`M_{gas}/M_{stars}=2.0`$, then becomes proportionately more disrupted with greater $`M_{gas}/M_{stars}`$ but survives (by our definition) to a $`M_{gas}/M_{stars}`$ of over $`4.0`$.

## 6 Discussion

We have modelled the gas expulsion from the core of a stellar cluster via a rapid change in the stellar velocities and thus an effective change in the gravitational potential of the core. This is taken to mimic the evolutionary stages of a hyper or ultra-compact HII region. Their evolution is likely to be driven by a stellar wind which stalls as mass loading occurs in order to explain their relatively long lifetimes $`(10^5yr)`$. Thus, we assume an instantaneous removal of gas from the core of the cluster but assume that the gas content of the rest of the cluster is for the most part unaffected during our simulation.

The formation of massive stars is largely an unsolved problem. What is certain is that they require high accretion rates in order to build the star up over timescales of $`10^5`$ to $`10^6`$ years. This is true whether they form just through accretion (Behrend & Maeder 2001; Bonnell & Bate 2002) or even through stellar mergers which require high accretion rates to reach the necessary stellar density (Bonnell et al. 1998; Bonnell & Bate 2002). The most obvious way that high accretion rates can occur is by situating the forming massive star in the centre of a stellar cluster. In this way, the overall cluster potential funnels matter down to the protomassive star where it can preferentially accrete it due to its position and mass (Bonnell et al. 2001). This naturally results in a mass segregated stellar cluster even while deeply embedded.

The simulations reported here show that gas expulsion is unlikely to remove this ’initial’ mass segregation unless the gas comprises more than 83 per cent of the total mass in the core. The core, being where the accretion rates are highest and the most massive stars are located, has therefore the lowest gas fraction of the system and is the most likely to be dominated by the mass in stars. Thus, the mass segregation which results from the cluster formation and gas accretion is likely to remain unaffected by the gas expulsion stage in the cluster.

In some young clusters such as the ONC, a few relatively massive stars are found to be some distance from the core of the cluster. If the above scenarios for the formation of massive stars are correct in that they form in the denser parts of the cluster, then we need to explain how they could have moved from the core to, for example, the half-mass radius of the cluster. the numerical simulations reported here offer one possibility. The massive stars could have originated in the core but escaped from this region during the gas expulsion phase (Williams 2002). We see that for gas fractions of greater than 83 per cent, some of the massive stars do escape the core and can be found virtually anywhere in the cluster. They will of course sink back to the cluster core but on longer timescales than the age of the cluster (Bonnell & Davies 1998).

The difficulty in the above is that relatively high gas fractions are required which may not be present in the core of the cluster. An alternative mechanism to explain the location of the massive stars not in the core would be stellar interactions with a binary or multiple system which eject one or more of the stars. Some of these ejections could be the high-velocity OB runaways while those that do not garner such a large velocity boost remain in the cluster but outside the core. If that is the case, their dynamics should resemble the massive stars ejected from the core in the simulations reported here.

## 7 Conclusion

Using a series of N-body simulations we follow the dynamical evolution of massive stars which are assumed to be initially located in the core of the cluster. The stars undergo evolution because of very rapid gas expulsion (assumed instantaneous) which is parametrised by the ratio of gas mass to stellar mass in the core $`(M_{gas}/M_{stars})`$. This parameter is varied from 0.25 to 5.0 (equivalent to a star formation efficiency of 80% down to 17%).

We find that the initial mass segregation in the core is destroyed for an initial gas mass of greater than 4 times the stellar mass, or a star formation efficiency of less than 20%. The initial conditions used are of a relaxed, mass-segregated cluster. The massive are therefore approximately in equipartition with low-mass stars in the core. Their lower velocities result in a lower probability of escaping the core once the gas too has been removed. This explains the relatively high value of $`M_{gas}/M_{stars}`$ required to dissolve the core. Otherwise, for star formation efficiencies of more than 20%, but less than 50% some massive stars will be ejected from the core and maybe beyond. However these ejected stars are found to have predominantly tangential motions, gaining angular momentum due to encounters along the way.


text:
# On susy periodicity defects of the Lamé Potentials

CONTRIBUTIONS REVISTA MEXICANA DE FÍSICA 49 SUPLEMENTO 2, 145–147 OCTUBRE 2003

## Abstract

The supersymmetric Quantum Mechanics approach is applied to embed bound states in the energy gaps of periodic potentials. The mechanism to generate periodicity defects in the first Lamé potential is analyzed. The related bound states are explicitly derived.

Let $`H=\frac{d^2}{dx^2}+V(x)`$ be a Schrödinger operator with $`V(x)`$ a real-valued and locally bounded function on $``$. The spectrum $`\sigma (H)`$ of $`H`$, as well as its resolvent set $`\rho (H)`$, will be then complementary subsets of $``$. The involved eigenvalue equation $`H\psi (x)=ϵ\psi (x)`$, admits solutions for arbitrary $`ϵ`$. The $`ϵ`$-values such that $`ϵ=E\sigma _d(H)`$ lead to square integrable wave-functions $`\psi _E(x)`$ in $`L^2()`$, while for $`ϵ\rho (H)`$ one gets eigenfunctions $`\psi (x)`$ deprived of physical meaning. It is well known that the class of Schrödinger’s operators $`H`$ admitting analytic expressions for the corresponding eigenfunctions is narrow. The supersymmetric (susy) Quantum Mechanics Coo00  intends to expand the class. In particular, the conventional Infeld-Hull factorization method has been extended to employ not only the physically interpretable $`\psi _E`$ but also the nonphysical eigenfunctions $`\psi `$ Mie84 . The main idea of this generalized factorization is to implement the Darboux-Crum transformation as the mechanism to construct multi-parametric families $`\{V(\lambda _1,\lambda _2,\mathrm{},\lambda _n,x)\}`$ of (almost) isospectral exactly solvable potentials everybody ; Fer97 . The members of each one of these families share common analytical properties, such as the number and type of (possible) singularities as well as the asymptotic behaviour. Hence, it is intended that the procedure does not add new singularities nor modifies the asymptotic behaviour of the initial potentials. This mechanism is ensured by taking eigenfunctions $`\psi `$ with $`ϵ\rho (H)`$ as the key transformation functions Fer97 .

The present contribution deals with the generalized factorization method applied to periodic potentials. As we shall see, the very interesting phenomenon of periodicity defects in crystal lattices can be produced as a consequence of the susy mechanism. We support our analysis mainly on the results reported recently in Fer02 .

Let us consider a Schrödinger operator $`H_0`$ for which the solutions $`u(x,ϵ)`$ to the eigenvalue equation

$$H_0u(x,ϵ)=ϵu(x,ϵ)$$
(1)
are known for any $`ϵ`$. Notice that equation (1) does not require the square integrability of $`u(x,ϵ)`$. An important example is the Hamiltonian with a real-valued, continuous and periodic potential

$$V_0(x+\tau )=V_0(x),x.$$
(2)
In this case, the translation operator $`Tu(x,ϵ)=u(x+\tau ,ϵ)`$, maps the two-dimensional subspace $`_ϵ`$ of solutions of the eigenvalue equation (1) into itself. As a consequence, there is a basis of eigenfunctions $`\{\mathrm{v}_1(x),\mathrm{v}_2(x)\}`$ in $`_ϵ`$ for which $`T`$ admits a $`2\times 2`$ matrix representation $`[T]`$ (the monodromy matrix). The eigenvalue equation for $`[T]`$, evaluated at $`x=T`$, determines two eigenvalues $`t_\pm `$ for a given energy $`ϵ`$, i.e., $`t_\pm =t_\pm (ϵ)`$. The corresponding eigenfunctions $`v_\pm (x,ϵ)`$ exist for any $`ϵ`$ and are called Bloch functions. The structure of $`v_\pm (x,ϵ)`$ depends essentially on the values of $`ϵ`$, namely, if $`ϵ`$ is such that $`t_\pm S^1`$, then $`ϵ\sigma (H_0)`$ and the related Bloch functions become complex and bounded. The spectrum of $`H_0`$ is then continuous, composed by a series of spectral bands in $``$. The edges of that bands are given by the $`ϵ`$-values such that $`t_\pm =\pm 1`$, and the Bloch functions become $`\tau `$-periodic or $`\tau `$-antiperiodic. Finally, if $`ϵ`$ is such that $`t_\pm /\{\pm 1\}`$, then $`ϵ\rho (H_0)`$ and the related Bloch functions are real and diverge at either $`\mathrm{}`$ or $`+\mathrm{}`$. Due to this mechanism, $`\rho (H_0)`$ splits into a (finite or infinite) sequence of forbidden gaps. In summary, $`\sigma (H_0)=\{ϵ;t_\pm S^1\}`$, while $`\rho (H_0)=/\sigma (H_0)`$.

As regards to the supersymmetric treatment of periodic Hamiltonians, one can consider two different options: (i) The Darboux transformations are such that, when applied to systems described by (2), the new potentials are also $`\tau `$-periodic. (ii) The procedure does not respect the global periodicity of the initial potential. In the former case, Dunne and Feinberg Dun98  have shown that, sometimes, the nonsingular first order transformed potential becomes simply a half-period displaced copy of the initial one. Such potentials were called self-isospectral. A further study (using the band edge eigenfunctions) produced the 2-susy self-isospectral potentials Fer00 . The above results have been recently generalized to embrace arbitrary $`\delta `$-displaced copies of the initial periodic potential Fer02 . Such a phenomenon has been named the translational Darboux invariance (or simply Darboux invariance), and it reproduces the self-isospectral case for $`\delta =\tau /2`$. The Darboux invariance is achieved by using Bloch solutions as transformation functions for $`ϵ`$ either at the band edges or in the resolvent set of the periodic Hamiltonian $`H_0`$. Concerning the option (ii), we first remark that arbitrary linear combinations of $`v_\pm (x,ϵ)`$ for $`ϵ\rho (H_0)`$, when used as transformation functions $`u(x,ϵ)`$ can produce regular non-periodic modifications in the periodic potential (2). As we are going to show, moreover, the functions $`u(x,ϵ)`$ allow to embed discrete energy levels in the gaps of $`H_0`$ and to determine the corresponding square integrable eigenfunctions Fer02 .

The generalized first order factorization method Mie84  allows to obtain the solutions $`\varphi _ϵ^{}(x,ϵ)`$ of a transformed eigenvalue equation $`H_1(ϵ)\varphi _ϵ^{}(x,ϵ)=ϵ^{}\varphi _ϵ^{}(x,ϵ)`$; $`H_1(ϵ):=\frac{d^2}{dx^2}+V_1(x,ϵ)`$, where the new potential $`V_1(x,ϵ)`$ is related with the initial one $`V_0(x)`$ by means of a Darboux transformation

$$V_1(x,ϵ)=V_0(x)+2\beta ^{}(x,ϵ).$$
(3)
The function $`\beta `$ here is nothing but the general solution of the Riccati equation

$$\beta ^{}(x,ϵ)+\beta ^2(x,ϵ)=V_0(x)ϵ.$$
(4)
At this stage, it is interesting to notice that the transformation $`\beta (x,ϵ)\frac{d}{dx}\mathrm{ln}u(x,ϵ)`$ allows to simplify the expressions for the potential (3) and the involved eigenfunctions $`\varphi _ϵ^{}(x,ϵ)`$ leading to:

$$\begin{array}{cc}V_1(x,ϵ)=V_0(x)2\frac{d^2}{dx^2}\mathrm{ln}u(x,ϵ),& \\ \varphi _ϵ^{}(x,ϵ)\frac{W(u,\psi _ϵ^{})}{u(x,ϵ)}& \end{array}$$

where, by simplicity, we have written $`uu(x,ϵ)`$ and $`\psi _ϵ^{}\psi _ϵ^{}(x)`$. The symbol $`W(u,\psi _ϵ^{})`$ represents the Wronskian of $`u`$ and $`\psi _ϵ^{}`$.

The analytical behaviour of $`\varphi _ϵ^{}(x,ϵ)`$ strongly depends on the values of $`ϵ`$ and $`ϵ^{}`$. If $`ϵ^{}=E\sigma (H_0)`$ and $`E_0>ϵ\rho (H_0)`$, with $`E_0\sigma (H_0)`$ denoting the lowest band edge, the eigenfunction $`\varphi _E(x,ϵ)`$ can be bounded. A particular detail of the method is that the function $`\varphi _ϵ(x,ϵ)1/u(x,ϵ)`$ satisfies identically $`H_1(ϵ)\varphi _ϵ(x,ϵ)=ϵ\varphi _ϵ(x,ϵ)`$ and becomes orthogonal to all the $`\varphi _E(x,ϵ)`$. Moreover, it can be chosen in $`L^2()`$. Thereby $`\sigma (H_1(ϵ))=\sigma (H_0)\{ϵ\}`$. In other words, the Hamiltonian $`H_1(ϵ)`$ can be constructed in order to have a bound state embedded in the lowest spectral gap of $`H_0`$. Hence, $`H_1(ϵ)`$ will be almost isospectral with $`H_0`$ for a properly selected $`ϵ`$. This situation allows to consider them as susy partner Hamiltonians. Concerning the analytical properties of $`V_1(x,ϵ)`$, it is very common to find that those conditions on $`ϵ`$, making $`\varphi _E(x,ϵ)`$ a square integrable function, lead to potentials (3) with exactly the same number and nature of divergences as the initial one (see for example Fer97 ). In particular, as for $`ϵ<E_0`$ the function $`u(x,ϵ)`$ can be chosen nodeless, the potential $`V_1(x,ϵ)`$ in (On susy periodicity defects of the Lamé Potentials) is just a regular non-periodic modification of $`V_0(x)`$. The Figure 1 shows one of these first order regular deformations for the simplest Lamé potential

$$V_0(x)=2m\mathrm{sn}^2(x|m).$$
(5)
Let us remark that the functions $`u(x,0.1)`$, producing the potential $`V_1(x,0.1)`$ in Figure 1, preserve the smoothness of (5) as well as its asymptotic behaviour. Moreover, the deformation of $`V_1`$ in the neighbourhood of $`x=0`$ indicates the presence of the bound state energy level at $`ϵ=0.1`$ in the first forbidden gap $`(\mathrm{},E_0)`$. For the present case, the related wave-function $`1/u(x,0.1)`$ is a smooth function of $`x`$, sharply peaked at $`x=0`$.

The procedure can be now iterated by considering $`ϵ=ϵ_a`$ as given and $`V_1(x,ϵ_a)`$ as the known exactly solvable potential. Thus, one has to look for a new potential $`V_2(x,ϵ_a,ϵ)`$ satisfying the second order relationship

$$V_2(x,ϵ_a,ϵ)=V_1(x,ϵ_a)+2\beta ^{}(x,ϵ_a,ϵ).$$
(6)
Therefore, the new function $`\beta (x,ϵ_a,ϵ)`$ must fulfill the Riccati equation $`\beta ^{}(x,ϵ_a,ϵ)+\beta ^2(x,ϵ_a,ϵ)=V_1(x,ϵ_a)ϵ`$, where $`ϵ`$ is arbitrary. The solutions of this equation are easily determined by knowing those of (4) for at least two different factorization constants $`ϵ`$, $`ϵ_a`$ Fer97 . Then, one can use the finite-difference Bäcklund algorithm Fer98 . It is remarkable that $`\beta (x,ϵ_a,ϵ)`$ appears to be symmetric under the interchange of $`ϵ_a`$ by $`ϵ`$ and vice versa. This means that one can achieve the potential (6) by finding first $`V_1(x,ϵ_a)`$ and then $`V_2(x,ϵ_a,ϵ)`$ or, first $`V_1(x,ϵ)`$ and then $`V_2(x,ϵ,ϵ_a)`$ Fer97 . More details concerning further iterations and symmetry properties of the Darboux transformations can be found in Nie01  and references quoted therein. Observe that potential (6) and the resulting eigenfunctions $`\varphi _ϵ^{}(x,ϵ_a,ϵ)`$ can be rewritten as:

$$\begin{array}{cc}V_2(x,ϵ_a,ϵ)=V_0(x)2\frac{d^2}{dx^2}\mathrm{ln}W(u_a,u);& \\ \varphi _ϵ^{}(x,ϵ_a,ϵ)\frac{W(u_a,u,\psi _ϵ^{})}{W(u_a,u)}& \end{array}$$

where, for the sake of simplicity, we have written $`u_au(x,ϵ_a)`$. Let us consider the precedent 1-susy step as given with $`ϵ_a<E_0\sigma (H_0)`$. The Hamiltonian $`H_1(ϵ_a)`$ is therefore regular and $`\varphi _ϵ^{}(x,ϵ_a,ϵ)`$ becomes an eigenfunction of $`H_2(ϵ_a,ϵ)`$ with eigenvalue $`ϵ^{}`$. If $`ϵ^{}=E\sigma (H_0)`$, $`ϵϵ_a`$ and $`E_0>ϵ\rho (H_0)`$, one easily finds that $`\varphi _E(x,ϵ_a,ϵ)`$ can be bounded. Yet, the function $`\varphi _ϵ(x,ϵ_a,ϵ)1/\varphi _ϵ(x,ϵ_a)`$ is also orthogonal to all the $`\varphi _E(x,ϵ_a,ϵ)`$ and can be chosen in $`L^2()`$. Hence $`\sigma (H_2(ϵ_a,ϵ))=\sigma (H_1(ϵ_a))\{ϵ\}`$ and the Hamiltonians $`H_1(ϵ_a)`$, $`H_2(ϵ_a,ϵ)`$ are almost isospectral. The iteration chain now establishes that $`H_2(ϵ_a,ϵ)`$ is also a susy partner of $`H_0`$ of the second order, i.e., $`H_2(ϵ_a,ϵ)`$ and $`H_0`$ are now “sususy partners”. In the present case, the iteration procedure leads to nonsingular Hamiltonians at each step, hence, the sususy transformation is said to be reducible.

Now, let us consider the case when $`E_0<ϵ_a\rho (H_0)`$. As the function $`u(x,ϵ_a)`$ in (On susy periodicity defects of the Lamé Potentials) has a finite number of zeros in $`(0,\tau )`$, the potential $`V_1(x,ϵ_a)`$, as well as the functions $`\varphi _ϵ^{}(x,ϵ_a)`$, become singular. However, there is still a chance to construct regular potentials $`V_2(x,ϵ_a,ϵ)`$ and functions $`\varphi _ϵ^{}(x,ϵ_a,ϵ)`$. It is enough to take $`ϵ_a,ϵ\rho (H_0)`$ such that $`E_j<ϵ_a,ϵ<E_j^{}`$, with $`E_j,E_j^{}`$ two subsequent band edges of $`H_0`$ (we have used the notation of Fer00 ). Hence, for an appropriate choice of the involved parameters, the Hamiltonian $`H_2(ϵ_a,ϵ)`$ is regular although the intermediary ones $`H_1(ϵ_a)`$ and $`H_1(ϵ)`$ are singular. We then say that the second order transformation is irreducible. Figure 2 shows the result of an irreducible second order regular transformation of potential (5) for $`m=0.5`$. The values of the $`ϵ`$’s are in the gap $`(E_1,E_1^{})=(1,1.5)`$. On the other hand, Figure 3 illustrates one of the related bound states for $`ϵ=1.45`$.

Finally, let us remind that the presence of bound states in the continuum is not a novelty in Quantum Mechanics. For example, Stillinger and Herrick used a modification of the von Neumann-Wigner approach to produce local potentials with bound eigenstates embedded in the continuum of scattering states Sti75 . Similar results are obtained by means of the susy Quantum Mechanics Pap93 . However, the situation when the bound states are embedded in the energy gaps of a periodic potential has been out of the attention for researches in the area. Further details on this subject can be found in Fer02 .

The support of CONACyT project 32086E and of CINVESTAV project JIRA’2001/17 is acknowledged.