format: text

text:
# The Fine Structure and Outskirts of DDO 154

## 1 Introduction

The nearby irregular galaxy DDO 154 has received a great deal of attention for its extraordinarily extended H I disk (Krumm & Burstein, 1984; Hoffman et al., 1993; Carignan & Purton, 1998) and its global dynamics, which seem to be dominated by dark matter at all radii (Carignan & Freeman, 1988; Carignan & Beaulieu, 1989) but require some physics beyond that of the Navarro, Frenk & White (1996, 1997) “universal” halo (Burkert & Silk, 1997; Gelato & Sommer-Larsen, 1999; Ghigna et al., 2000; Milgrom & Braun, 1988). Since the H I disk extends so far beyond the visible stars while the galaxy is apparently isolated from all neighbors, in contrast to most other examples of H I envelopes that extend far outside their optical components, DDO 154 would seem to be an excellent case in which to study the structure and thermodynamic state of a relatively pristine gas cloud. There are two major questions: (1) How does the neutral gas disk terminate — is it truncated abruptly by the extragalactic UV radiation field (Corbelli & Salpeter, 1993; Maloney, 1993; Dove & Shull, 1994), by UV photons leaking out from the central stellar disk (Bland-Hawthorn, Freeman & Quinn, 1997), or by some other process? (2) What is the fine structure and velocity dispersion of the outer gas disk — is that gas, presumably well removed from energy input from star formation and supernovae, smoothly distributed or as filamentary as the inner parts of the Magellanic Clouds (Stanimirović et al., 1999; Kim et al., 1998) or well resolved spiral galaxies (Braun, 1997)?

To broach (1), we extended our Arecibo map of the outskirts of the galaxy, taking full advantage of the enhanced sensitivity and greatly reduced sidelobes of the upgraded telescope. For (2), we undertook to improve, by a factor of three, the angular resolution of the mapping of the H I envelope to assess the fine structure of the gas disk and in the hopes of finding sufficient column density to perform emission/absorption studies against one or more of the NVSS continuum sources in the outskirts of the SW quadrant of the H I disk. The average column density from the D array, DRAO and Arecibo mapping cited above would be too low, but if the outer envelope proved to be filamentary there would be a reasonable probability that such studies would be feasible.

In Sect. 2 our Arecibo study of the outskirts of the southernmost quadrant of the galaxy is presented. We detail our VLA C array map of the galaxy in Sect. 3. Sect. 4 is a discussion of the various theoretical points of interest, and a summary and conclusions follow in Sect. 5.

## 2 Arecibo Observations

Our mapping of the outskirts of the southern quadrant of DDO 154 was conducted as a commissioning phase project at the upgraded Arecibo Observatory 111The Arecibo Observatory is part of the National Astronomy and Ionosphere Center, which is operated by Cornell University under a management agreement with the National Science Foundation. in July/August 1998. We used the Gregorian feed system with the “L narrow” receiver in total power (position-switched) mode, with 6.10 kHz (about 1.3 km $`\mathrm{s}^1`$) channel spacing. Calibration was accomplished by observing several continuum sources from the VLA calibrator list, chosen to have small size compared to the $`3\stackrel{}{\mathrm{.}}2`$ beam. In addition, we reobserved several spiral galaxies for which we had high signal-to-noise pre-upgrade H I measurements and which were known to be $`3\stackrel{}{\mathrm{.}}2`$ in extent.

At the time of these observations, final focussing of the Gregorian feed system was still underway and the first sidelobe ring at 21cm was less uniform than in the final configuration and varied significantly with azimuth and zenith angle. For a fuller characterization of the first sidelobe ring of the current instrument, see Heiles et al. (2001). By scanning the beam across strong point sources at close to DDO 154’s declination as they drifted across the field of view, we were able to determine that we could keep the center of the galaxy in the weakest part of the sidelobe ring (15dB or more down from the peak of the beam) as long as the point being observed lay between the center of the galaxy and the center of the dish at that moment. Consequently, we planned our observations day by day to make sure that we kept that orientation, starting with points on the SE side of the galaxy as it entered the field of view and working around the rim to points on the SW side as the galaxy left the field of view.

The sidelobe contributions are not negligible, however, even though they are greatly reduced from those of the pre-upgrade circular feed. The beam had a null ring (at least 20dB down) about $`9\mathrm{}\times 7\stackrel{}{\mathrm{.}}5`$ in diameter, and the first sidelobe peaked in a ring about $`11\stackrel{}{\mathrm{.}}6\times 10\stackrel{}{\mathrm{.}}6`$ in diameter with a maximum that varied from about $`16`$dB to about $`13`$dB; as discussed above, only the weaker part of the ring extended toward the center of the galaxy in our observations. A $`15`$dB uniform ring would have a total effective area 25-30% of the main beam, so that one sextant of the ring would give about 4% of the response of the main beam. To correct for those sidelobe contributions, we interpolated between the points on our pre-upgrade flat feed map of the inner portions of the galaxy (Hoffman et al., 1993) to estimate the flux that would fall in each sextant of the sidelobe ring, and then subtracted 4% of that flux from the spectrum obtained at the position of the main beam. The spectrum obtained at each observed position and the estimated sidelobe contribution to that spectrum are displayed in Figure 1. The positions in relation to those mapped with the pre-upgrade flat feed in Hoffman et al. (1993) are shown in Fig. 2.

There is an additional worry: scattering from the triangular platform and support cables is expected to produce far sidelobes of uncertain magnitude (but considerably less than the first sidelobe ring) which vary greatly with azimuth and zenith angle. That expectation is consistent with our observation that the ability of the first sidelobe ring to account for the observed flux varies erratically from point to point, at least away from the extension of the major axis. Far sidelobes seem to be a more reasonable explanation for that patchy low-level H I, since the emission not accounted for by the first sidelobe has a velocity width comparable to that of spectra from points near the galaxy’s center, not as narrow as the outer VLA C array spectra shown in the following section. For the far sidelobes to account for the residual flux we see, the magnitude of the response for the portion of the far sidelobe which catches the center of the galaxy would have to be about $`25`$dB for a wide range of distances from the center of the beam. Nevertheless, in those cases where the remnant flux falls within the velocity bounds of the profiles obtained near the center of the galaxy, we conservatively attribute that flux entirely to incompletely removed sidelobes. The only points with marginally significant emission that we cannot so attribute to sidelobes are those at the end of the warped major axis, labeled “J” and “K,” since that flux is seen at velocities higher than any from the central parts of the galaxy. Taking, conservatively, the uncertainty in the corrected flux to be $`1/3\times 2\times (50\mathrm{km}\mathrm{s}^1)\times (\mathrm{rms}\mathrm{flux})`$ added in quadrature to the uncertainty in the sidelobe contribution (approximately the square root of the number of contributing points times 20% of the sidelobe flux), the sidelobe-corrected emission we detect at points J and K is marginally significant at levels of $`2\sigma `$ and $`3\sigma `$, respectively. In Fig. 3 we show the sum of those two spectra. Calculating uncertainties as above, the summed signal is significant at just over $`3\sigma `$.

From the measured fluxes, after subtraction of the first sidelobe contributions, H I column densities can be calculated. The profiles obtained for points in this quadrant of the galaxy are approximately $`50\mathrm{km}\mathrm{s}^1`$ broad, so in those cases where the sidelobe subtraction has left $`<2\times (50\mathrm{km}\mathrm{s}^1)\times (\mathrm{rms}\mathrm{flux})`$ we take that expression to be the effective upper limit on the corrected flux. The fluxes and column densities, both as measured and corrected, are given in Table 1.

The (warped) major axis of the galaxy follows points 1-4, J and K through the SW quadrant. A second radial spoke, for comparison, can be taken to include points y, z, E and F. The galaxy’s center is at the point of convergence of the two spokes, marked by a filled circle. The run of column density (corrected for first sidelobe contamination) vs. radius along those two spokes is shown in Fig. 4. For the major axis, the center and points 1 through 4 are from our pre-upgrade flat-feed data as reported in Hoffman et al. (1993). Uncertainties in these points are of order 10%. The outer two points are from our new data (J and K) for the SW end of the major axis. For these, at each radius we display two points connected by a vertical bar. The upper point ignores sidelobe contamination and should be regarded as an absolute upper limit to the column density at that radius. The lower point is corrected for sidelobe contamination as discussed above. On these two particular spokes, the sidelobe-corrected fluxes exceed $`2\times (50\mathrm{km}\mathrm{s}^1)\times (\mathrm{rms}\mathrm{flux})`$ at every point. For the off-major axis, however, we cannot rule out incompletely removed sidelobes as the source of the flux as discussed above. Therefore we are led to interpret these results as indicating a truncated edge to the neutral hydrogen disk at a radius of about 15′ and a column density around $`10^{19}\mathrm{atoms}\mathrm{cm}^2`$, with marginal evidence for the column density levelling out at $`2\times 10^{18}\mathrm{atoms}\mathrm{cm}^2`$ beyond that point on the major axis. The uncertainty in the outermost column density is on the order of a factor of 2, however, and the conservative viewpoint would be that we have only an upper limit on the column density at point K.

## 3 VLA C Array Mapping

H I spectral line mapping using 27 antennas in the C array of the Very Large Array 222The Very Large Array of the National Radio Astronomy Observatory is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc. was conducted on 24 Dec 1998. Observational details are given in Table 2. The pointing center was about 7′ SW of the center of the galaxy, to give us greater sensitivity in the region where the H I is most extended and where the NVSS sources are concentrated. To obtain information about the velocity dispersion in the outer parts of the H I envelope, we used 2.5 km $`\mathrm{s}^1`$ velocity channels. Online Hanning smoothing was employed, and calibration was accomplished using sources 1219+285 (B1950) and 3C286 from the VLA calibrator list. The data were calibrated and editted using standard tasks in the Astronomical Image Processing System (Classic AIPS). Continuum subtraction was done in the $`uv`$ plane via UVBAS, and maps were made and CLEANed using IMAGR with zero-spacing fluxes estimated from our Arecibo map. Two cubes were made, one with robustness set equal to 0 (beam $`18\mathrm{}\times 13\mathrm{}`$) and one with robustness 1 (beam $`18\mathrm{}\times 17\mathrm{}`$). After imaging, each data cube was corrected for the VLA primary beam. For quantative results we have used the robustness 1 cube, since it is slightly more sensitive. For moment analyses we blanked the cube using a mask made by smoothing the cube to a $`60\mathrm{}\times 60\mathrm{}`$ beam, blanking everything more than 20 pixels outside the mask in each plane.

We detect H I out to a radius of about 6′ from the center of the galaxy, to a column density detection limit of $`4\times 10^{19}\mathrm{atoms}\mathrm{cm}^2`$ when the unblanked cube is smoothed to the typical linewidth ($`19\mathrm{km}\mathrm{s}^1`$) of features within individual beams in the outer part of the H I disk. We conclude that the outermost gas is quite smooth on this scale, with no clumps that exceed $`1\times 10^{20}\mathrm{cm}^2`$ on the scale of our $`18\mathrm{}\times 17\mathrm{}`$ beam, or $`280\times 260`$ pc assuming a 3.2 Mpc distance to the galaxy (Carignan & Beaulieu, 1989). The more diffuse emission detected in our Arecibo map (Hoffman et al., 1993) and in the D array plus DRAO map of Carignan & Purton (1998) has been resolved out by the lack of short baselines, and of course we can say nothing about structure on scales smaller than our beam size. There are no indications of any region $`18\mathrm{}`$ or more in size with sufficient column density overlying any of the NVSS sources for absorption studies, unfortunately. However, there is considerable structure in the H I as far out as we can detect gas at this resolution, as shown in Fig. 5. There the total H I emission, integrated over velocity after blanking as described above, is shown. It does not reach to our minimum detectable column density since blanking cannot be done perfectly; however, we have inspected a cube made by integrating over sets of 5 successive channels without blanking and found no features significantly outside the emission shown in Fig. 5.

### 3.1 Two Prominent Cavities

In particular, there are two prominent cavities centered approximately at $`12^h51^m33^s+27\mathrm{°}24\mathrm{}20\mathrm{}`$(1950) and $`12^h51^m44^s+27\mathrm{°}24\mathrm{}50\mathrm{}`$(1950). In Fig. 6 each channel map is shown with a single contour indicating the outermost extent of the optical images from the Digitized Sky Survey (DSS). The more easterly cavity lies on the very edge of the optical image, just beyond the 26.5 mag $`\mathrm{arcsec}^2`$ $`B`$ isophote as shown in Carignan & Beaulieu (1989). In the DSS image, there is a faint smudge that might be an H II region near the center of the cavity. The more westerly cavity, however, has no indication of any starlight whatsoever within it. Neither can any H$`\alpha `$ emission be found associated with either cavity in the image acquired by Hunter, Elmegreen & Baker (1998), reproduced in Fig. 7, nor in the deeper H$`\alpha `$ image of Kennicutt & Skillman (2001).

The diameters of the cavities are about 40″and 60″, respectively, or about 600 pc and 900 pc for our assumed 3.2 Mpc distance to DDO 154. The velocity depths are about 20.6 and 33.6 km $`\mathrm{s}^1`$ (corresponding to expansion velocities of order 10.3 and 16.8 km $`\mathrm{s}^1`$). However, as shown in Figs. 8 and 9 the spectra in the interiors of the cavities do not give clear indications of expanding walls at the front and rear; rather, the gas in the interior seems to partake of the same velocity field as the gas immediately adjacent to the cavity — the neutral component in each cavity is simply less dense. There are no indications of vorticity, although in the major axis-velocity map (Fig. 10) there is a dense ridge at the inner edge of the western cavity which clearly departs from uniform circular motion.

The velocity field contours from fitting gaussians to our C array mapping are shown in Fig. 5. A position-velocity map, along the major axis and summed over the minor axis, is shown in Fig. 10. The more easterly and more westerly holes fall at about 0 and $`80\mathrm{}`$, respectively, on the position axis, and there is a dense ridge of material that deviates from circular motion just inside the westerly hole. The rotation curve derived from the velocity field is shown in Fig. 11. It does not differ significantly in the inner parts from the lower resolution mapping of Carignan & Freeman (1988) or Carignan & Purton (1998), so evidently beam-smearing (Swaters, Madore & Terwhella, 2000; van den Bosch et al., 2000) has not significantly affected the discussions of mass modeling cited above. Further analysis and discussion of the rotation curve will be deferred to a forthcoming paper (Swaters et al., 2001).

The enhanced resolution in angle and velocity of our observations do allow us to determine the velocity dispersion of the gas. Selected velocity profiles are displayed in Fig. 12; each is fit well by a single gaussian with no significant residuals. A map of FWHM of gaussians fit to profiles at each pixel is presented in Fig. 13. The corresponding dispersions are 7 to 8 $`\mathrm{km}\mathrm{s}^1`$ over most of the outer parts of the galaxy, rising to about 9 $`\mathrm{km}\mathrm{s}^1`$ in the central part (more or less coincident with the part that has formed stars). The rising portion of the rotation curve would increase the velocity dispersion in the central regions (out to radius $`50\mathrm{}`$) of the galaxy by somewhat less than 2 km $`\mathrm{s}^1`$, about what is observed, so it is not clear whether any part of the increase in the dispersion in the center can be attributed to heating by stars. What is clear is that the dispersion in the outer parts, well outside the stellar portion, is still characteristic of the warm neutral medium (gas temperature of order $`10^4`$ K) rather than anything much colder.

The scaleheight of the gas in the outer parts of the disk can be estimated following Puche et al. (1992). We adopt $`7.5\mathrm{km}\mathrm{s}^1`$ as a typical velocity dispersion at radii $`4`$$`6\mathrm{}`$. The volume mass density in the disk at that radius is essentially all due to dark matter (Carignan & Purton, 1998), which we take to be distributed spherically for convenience. Then the mass density at radius $`r`$ is just

$$\rho (r)=\frac{1}{4\pi Gr^2}\frac{d}{dr}\left[v(r)^2r\right]$$
(1)
where $`v(r)`$ is the circular velocity at radius $`r`$. Where the rotation curve is flat, this reduces to

$$\rho (r)=\frac{v(r)^2}{4\pi Gr^2}.$$
(2)
In the case of DDO 154, the rotation curve goes flat at a radius of $`4.5`$ kpc (Carignan & Purton, 1998), then declines gradually. Taking the circular velocity to be $`47\mathrm{km}\mathrm{s}^1`$ at that point, we find a volume density $`\rho =1.37\times 10^{25}\mathrm{g}\mathrm{cm}^3=0.002\mathrm{M}_{\mathrm{}}\mathrm{pc}^3`$ at that radius. The scaleheight is

$$h=\frac{\sigma }{[4\pi G\rho ]^{1/2}}$$
(3)
(Kellman, 1972; van der Kruit, 1981), giving $`h=718`$ pc for DDO 154 at a radius of 4.5 kpc ($`4.\mathrm{}8`$), well outside the stellar disk.

## 4 Discussion

### 4.1 The outer envelope

The most reasonable interpretation of the Arecibo results presented in Sect. 2 is that the H I phase truncates rather abruptly, falling by nearly an order of magnitude in less than $`3\stackrel{}{\mathrm{.}}5`$ (3.2 kpc) at a column density around $`10^{19}\mathrm{atoms}\mathrm{cm}^2`$, with at best marginal indications of emission an order of magnitude lower continuing along the major axis. As discussed in Sect. 2, the uncertainty in the column density at point J is of order a factor of 2, and the outermost point (K) should perhaps be regarded as an upper limit only. While the suggested upturn in the rotation curve at large radii is intriguing, the data in hand do not warrant much further discussion of the marginally significant continuation of the major axis. Still more Arecibo mapping, and a detailed study of the sidelobes of the upgraded Arecibo beam, is indicated.

That the neutral phase should abruptly decrease at such a column density is no surprise, since similar truncation has been observed in M33 (Corbelli, Schneider & Salpeter, 1989; van Gorkom, 1993). A theoretical explanation, in terms of photoionization by extragalactic UV, has been offered by Corbelli & Salpeter (1993); Maloney (1993); Dove & Shull (1994). However, detection of H$`\alpha `$ flux around NGC 253 by Bland-Hawthorn, Freeman & Quinn (1997) is claimed by those authors to require a more local photoionizing source, e.g., UV leaking out of the disk of the galaxy. DDO 154, forming stars much more quiescently, is unlikely to be leaking enough H$`\alpha `$ to produce the truncation we claim to see so far outside the star-forming region. A definitive test would be to compare the truncation along the line of nodes of the warped H I disk with the direction perpendicular to that, since UV would be expected to escape mainly perpendicular to the disk of the galaxy and so could truncate the disk where a warp has lifted the edge of the disk out of the central plane but not along the line of nodes (Bland-Hawthorn, Freeman & Quinn, 1997). Since we see, if anything, more abrupt truncation along the off-major axis spoke, which should be closer to the line of nodes, the data presented here suggests that the photoionizing flux is extragalactic in this case.

From our Arecibo data (Hoffman et al., 1993, and Sect. 2) and the VLA + DRAO map of Carignan & Purton (1998) we know that H I remains above a few $`\times 10^{19}\mathrm{atoms}\mathrm{cm}^2`$ for another 6 or 8′ beyond that detected in the data of Sect. 3. The outer gas must therefore be rather smooth, without structure that would amount to density contrasts of more than a factor of 2 or 3 on the scale of our 18″ beam, since our C array detection threshold is a bit less than $`10^{20}\mathrm{atoms}\mathrm{cm}^2`$ on that scale and much emission at $`\mathrm{few}\times 10^{19}\mathrm{atoms}\mathrm{cm}^2`$ is evidently resolved out. Structure such as we see in the inner 6′ , with density contrasts of a factor of 10 on 18″ scales, does not continue into the outer envelope. On the surface, this argues for the structure ultimately to be generated by stellar processes. However, we argue in the following subsection that the connection between stars and H I structure cannot be too intimate.

### 4.2 Cavity generating mechanisms

The cavities described in Sect. 3.2 are fairly typical of the holes found in Holmberg II (Puche et al., 1992), IC 2574 (Walter & Brinks, 1999), and DDO 47 (Walter & Brinks, 2001), among others. A similar cavity can be seen in the channel maps for NGC 2366 (Hunter, Elmegreen & van Woerden, 2001). The required energies are consistent with those provided by multiple supernovae in a star formation region. However, at least one of these cavities is well outside the region where any star formation has taken place. Nothing is visible in the interior of the western cavity in B, R or H$`\alpha `$ images. Similar claims have been made for many of the H I holes in Holmberg II (Rhode et al., 1999), although those holes do not appear to be as far removed from the stellar regions of that galaxy. Imaging in the UV might be more definitive (Stewart & Walter, 2000), but to our knowledge no UV image with sufficient resolution has been published. There is considerable resemblance to the giant H I hole in NGC 6822 (de Blok & Walter, 2000), although DDO 154 gives no indication of tidal interaction such as that claimed for NGC 6822.

Another hole-making mechanism that has been discussed (Santillán et al., 1999, and references therein) is a collision with a High Velocity Cloud. We cannot rule that out in this case, although we otherwise have no evidence of there being HVC around DDO 154 and would not expect high impact velocities since the galaxy has low mass. A third mechanism has been proposed recently by Crosthwaite, Turner & Ho (2000)) and supported by numerical hydrodynamical simulations (Wada, Spaans & Kim, 2000): cavities formed as gravitational instabilities in the process of spiral arm formation. While there are no clear indications of spiral arms in either the optical or H I images of DDO 154, a gravitational instability of this sort may still be the best explanation of these cavities. That would require, however, that the bulk of the mass of the galaxy be in the disk, distributed spatially as the H I is distributed. Other researchers have found that the best fit to the VLA + DRAO rotation curve has dark matter in a slightly flattened ($`b/a=0.6`$) halo (Carignan & Purton, 1998) or in a disk as flat as the H I disk (Hoekstra, van Albada & Sancisi, 2001), in each case with density proportional to that of the H I. Modified Newtonian Dynamics (MOND) has been claimed to give a good fit to the rotation curve (Milgrom & Braun, 1988) with no dark matter at all.

## 5 Conclusions and summary

To study physical processes in the outermost H I of a galaxy that is relatively undisturbed by tidal influences, we have pursued complementary mapping of DDO 154 with the Gregorian feed system at Arecibo Observatory and with the C array of the Very Large Array. The $`3\stackrel{}{\mathrm{.}}2`$ Arecibo beam gives us the most sensitive available information about the outermost edges of the H I disk, and the C array mapping gives us information about the fine structure of the gas at somewhat smaller radii, but still well outside the stellar disk.

Despite uncertainty about the contribution of far sidelobes to our Arecibo spectra, we interpret our results to show a truncation of the H I disk of DDO 154 at a column density around $`10^{19}`$ atoms $`\mathrm{cm}^2`$, consistent with the truncations observed around the spirals M33 (Corbelli, Schneider & Salpeter, 1989) and NGC 3198 (van Gorkom, 1993). This is consistent with theoretical expectations (Corbelli & Salpeter, 1993; Maloney, 1993; Dove & Shull, 1994) for the truncation produced by the extragalactic UV field. There are tantalizing hints that the column density along the continuation of the major axis might level off near $`2\times 10^{18}`$ atoms $`\mathrm{cm}^2`$ and at higher velocity than is seen in any other spectrum for the galaxy, but with only marginal significance. A detailed study of the sidelobes of the Gregorian feed system at Arecibo, and further mapping of DDO 154 with the Arecibo beam, will be required to verify these results.

The VLA results show that the gas beyond $`6\mathrm{}`$ in radius must be relatively smooth on scales of a few hundred pc or larger, with no structure at those scales exhibiting a density contrast of a factor of 10 or more. We do not find sufficient column density against any background continuum sources to warrant absorption studies. However, there is considerable structure in the gas disk at smaller radii, even well outside the regions where there are visible stars. Two prominent cavities well removed from any significant stellar populations were noted, with a pronounced overdensity between them. While the energies required for evacuation are consistent with those produced by multiple supernovae, there is no visible trace of the stars which could have engendered those supernovae within a kpc of the center of the larger cavity. The smaller of the two cavities is centered just outside the 26.5 mag $`\mathrm{arcsec}^2`$ $`B`$ isophote. Other hole-making mechanisms were discussed in Sect. 4.2, but none presents itself as the obvious choice.

The velocity dispersion of the gas, measured within our $`18\mathrm{}`$ (270 pc) beam, is 7 to 8 km $`\mathrm{s}^1`$ throughout the disk (to $`6\mathrm{}`$ radius), with a slight enhancement at the galaxy’s center (which could be produced either by increased turbulence associated with star formation or by the rising part of the rotation curve) and a slight diminishment in the center of the larger cavity. Velocity profiles are well fit by single gaussians at all points. This translates to a disk scaleheight of 718 pc at the radius where the rotation curve turns flat, about 4.5 kpc.

Detailed discussion of the rotation curve will follow in a separate paper (Swaters et al., 2001).

D. Hunter kindly made available to us her H$`\alpha `$ image of DDO 154. We thank J. van Gorkom for guidance on observing procedures and suggestions on analysis, and R. Swaters and R. Braun for valuable discussions. Partial support was provided by the Crafoord Fund at Cornell University and by a faculty research grant from Lafayette College. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope, which is operated by the California Institute of Technology and Palomar Observatory on Palomar Mountain. The plates were processed into the present compressed digital form with the permission of that institution. The Second Palomar Observatory Sky Survey (POSS-II) was made by the California Institute of Technology with funds from the National Science Foundation, the National Geographic Society, the Sloan Foundation, the Samuel Oschin Foundation, and the Eastman Kodak Corporation. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.


text:
# Non-sequential triple ionization in strong fields

## I Introduction

Multiple ionization is a fundamental process in laser-atom interactions . It received particular attention when it was realized  that the cross section was much larger than could be expected for independent electrons model . One, therefore, came to distinguish two modes of multi-ionization, a sequential ionization processes that is compatible with an independent electron model, and a non-sequential one, where correlations in electron dynamics are important. The way interactions come in is a matter of debate but the rescattering model  has a fascinating combination of quantum tunneling and classical dynamics of electrons in a field and has the strongest support from experimental observations, theoretical analyses and numerical simulations . One consequence of the rescattering model is that the final decay towards multiple ionization is preceeded by the formation of a highly excited complex of electrons close to the nucleus. This compound states forms during the collision of the rescattered electron with the core. In our analysis we take this state as the initial condition for the final decay towards a multiply ionized atom.

The way two electrons escape from this compound state in non-sequential double ionization has been elucidated in a recent set of experiments by Weber et al . They found that the electrons escape preferentially with the same momenta. This observation prompted us to analyze the classical pathways that could lead to double ionization . In an extension of Wanniers analysis  we found that for two electron escape most classical paths lie near a subspace of electrons moving with the same distance from the nucleus but reflection symmetric with respect to the field axis. Calculations of ion momenta parallel and perpendicular to the field within this subspace  are in favorable agreement with experiment.

At these laser intensities also more than two electrons can escape simultaneously. In particular, triple ionization has been investigated by Moshammer et al.  and the recoil ion momentum distributions have been measured. It is our aim here to present the extension of the previous arguments  to the case of triple ionization. We identify saddles for non-sequential three electron escape and the relevant symmetry subspaces in which the saddles are situated. Numerical simulations within the subspaces allow us to obtain the ion momenta that result from simultaneous triple ionization. Specifically, in section II we discuss the pathways that lead to triple ionization, identify two saddle configurations and analyze their stability properties. In section III we then present numerical results for the distributions of ion momenta within these two subspaces and compare to the experimentally observed distributions. Finally, in section IV we draw conclusions and given an outlook to higher multiple ionization.

## II Pathways to non-sequential triple ionization

As in the case of double ionization we can assume that the electron dynamics during the final decay is fast compared to the changes in the field . Then the straightforward extension of the double ionization analysis  to three electron escape calls for an investigation of stationary configurations of three electrons in an atom exposed to a strong static field. Keeping a high symmetry between the electrons an obvious candidate is a configuration where electrons are placed at the vertices of an equilateral triangle whose plane is perpendicular to the field axis (see Fig. 1). Deviations from this plane, in particular displacements in field direction, will again be amplified by electron repulsion, thus leading to electrons being pushed towards the nucleus and not to triple ionization. This configuration has $`C_{3v}`$ symmetry and will be analyzed in subsection II A.

Further investigations show, however, that there is a second stationary configuration. This second configuration has all electrons in a plane with one electron on the field axis and the others symmetric with respect to it. This configuration is of symmetry $`C_{2v}`$ and will be analyzed in subsection II B. The decay signatures of both configurations are different and the experimental observations can be used to decide between them.

In both cases we start from the Hamiltonian for three electrons in a neutral atom in a laser field polarized along the $`z`$-axis. Within dipole approximation, with infinitely heavy nucleus and in atomic units, the Hamiltonian is

$$H=\frac{1}{2}(𝐩_1^2+𝐩_2^2+𝐩_3^2)+V(𝐫_1,𝐫_2,𝐫_3,t),$$
(1)
with the potential energy

$`V`$ $`=`$ $`{\displaystyle \frac{3}{r_1}}{\displaystyle \frac{3}{r_2}}{\displaystyle \frac{3}{r_3}}`$ (4)
$`+{\displaystyle \frac{1}{|𝐫_1𝐫_2|}}+{\displaystyle \frac{1}{|𝐫_2𝐫_3|}}+{\displaystyle \frac{1}{|𝐫_1𝐫_3|}}`$
$`(z_1+z_2+z_3)Ff(t)`$
and the time dependence of the pulse

$$f(t)=\mathrm{sin}^2(\pi t/T_d)\mathrm{cos}(\omega t+\varphi ).$$
(5)
$`F`$, $`T_d`$, $`\omega `$ and $`\varphi `$ stand for peak amplitude, duration, frequency and phase of the external field, respectively.

The projection of total angular momentum onto the polarization axis is conserved so that the number of degrees of freedom is reduced by one, leaving eight plus the time dependence from the external field. We work in a subspace where this component vanishes.

Before we proceed with our analysis note a scaling symmetry of the classical Hamiltonian (1) that can be used to eliminate one parameter. If the variables are rescaled according to

$`H`$ $`=`$ $`F^{1/2}H^{}`$ (6)
$`𝐫`$ $`=`$ $`F^{1/2}𝐫^{}`$ (7)
$`𝐩`$ $`=`$ $`F^{1/4}𝐩^{}`$ (8)
$`t`$ $`=`$ $`F^{3/4}t^{}`$ (9)
$`\omega `$ $`=`$ $`F^{3/4}\omega ^{}`$ (10)
the dynamics becomes independent of the peak field amplitude, i.e. the system is described by the Hamiltonian (1) with $`F=1`$. We will use this scaling, but drop the prime to denote the scaled variables, for the analysis of the saddles, but we will keep the full field dependence for the dynamical simulations.

### A $`C_{3v}`$ symmetry subspace

For zero total angular momentum projection on the polarization axis, the $`C_{3v}`$ symmetric configuration corresponds in cylindrical coordinates to electron positions given by $`z_i=Z`$, $`\rho _i=R`$ and $`\phi _i=2\pi i/3`$ with conjugate momenta $`p_{z_i}=p_Z/3`$, $`p_{\rho _i}=p_R/3`$ and $`p_{\phi _i}=0`$, respectively. The Hamiltonian of the system in the $`C_{3v}`$ subspace then reads

$$H(p_R,p_Z,R,Z,t)=\frac{p_R^2+p_Z^2}{6}+V(R,Z,t),$$
(11)
with potential energy

$$V=\frac{9}{\sqrt{R^2+Z^2}}+\frac{3}{2R\mathrm{sin}(\pi /3)}3Zf(t)$$
(12)
The three terms in Eq. (12) are the cumulative interaction of three electrons with the triply-charged nucleus, the repulsion energy between electrons at distances $`2R\mathrm{sin}(\pi /3)`$ and the interaction with the external field, respectively. The variables are scaled according to (10) so that the field amplitude $`F`$ is absent.

Electron motion close to the nucleus is much faster than the change of the phase of the laser field applied in the experiments . Therefore, we can use an adiabatic approximation and gain insight into the qualitative features of the ionization process by analyzing the potential (12) for fixed external field. Note, however, that we use the full time dependence for the determination of the final ion momenta below. The potential (12) for a given time and thus a fixed value of $`f(t)`$ is shown in Fig. 2. The saddle is located along the lines $`Z_s=r_s\mathrm{cos}\theta _s`$ and $`R_s=r_s\mathrm{sin}\theta _s`$ with $`\theta _s=\theta `$ or $`\theta _s=\pi \theta `$ where

$$\theta =\mathrm{arctan}\frac{1}{\sqrt{2}}35^{}$$
(13)
and

$$r_s^2=\sqrt{6}/|f(t)|.$$
(14)
The energy of the saddle is

$$V_s=6^{3/4}2\sqrt{|f(t)|}.$$
(15)
During a field cycle the saddle moves in from infinity along the line $`\theta _s=\theta `$, back out to infinity and then in and out again along the line $`\theta _s=\pi \theta `$.

For the dynamics in full phase space the stability properties of the configuration are important. Since the electrons move fast compared to the frequency of the field we can for this analysis employ an adiabatic approximation and work with a constant field. Then we can scale out the field amplitude completely and put $`f(t)=1`$. The triangular electron configuration then becomes stationary and the dynamics near this saddle is governed by the eigenvalues of the linearization. Within the $`C_{3v}`$ symmetry subspace the triangular electron configuration has one unstable direction corresponding to a crossing of the saddle. In these scaled units the Lyapunov exponent for motion along this reaction coordinate  is $`\lambda _r=1.1054`$. The other eigenmode of the saddle in the subspace is stable.

Harmonic approximation in the full eight-dimensional configuration space results in three additional pairs of degenerate eigenvalues. Of these pairs two are stable and one is unstable. The unstable eigenspace points towards asymmetric configurations, has a Lyapunov exponent $`\nu _a=1.4496`$ and is spanned by the vectors

$`\rho _1=R_s+w_a,`$ $`z_1=Z_s+2.0642w_a,`$ (16)
$`\rho _2=R_s+w_a,`$ $`z_2=Z_s+2.0642w_a,`$ (17)
$`\rho _3=R_s2.0000w_a,`$ $`z_3=Z_s4.1284w_a,`$ (18)
and

$`\rho _1=R_su_a,`$ $`z_1=Z_s2.0642u_a,`$ (19)
$`\rho _2=R_s+u_a,`$ $`z_2=Z_s+2.0642u_a,`$ (20)
$`\rho _3=R_s,`$ $`z_3=Z_s.`$ (21)
The first mode describes a motion where for positive $`w_a`$ two electrons move away from and one towards the nucleus or vice versa for negative $`w_a`$. The other mode leaves one electron on the saddle and brings one closer and one further out.

One may compare the period of the laser field applied in the experiments  with the time scales for crossing of the saddle in the $`C_{3v}`$ subspace and for departure from the subspace. For peak field intensity $`1.510^{15}`$ W/cm2 and for wavelength of 795 nm, the field period (in the scaled variables (10)) is $`2\pi /\omega =33.6`$ while $`1/\lambda _r=0.9`$ and $`1/\nu _a=0.69`$. This indicates that crossing the saddle and departure from the non-sequential triple ionization manifold take place rather quickly, justifying the adiabatic analysis a posteriori.

### B $`C_{2v}`$ subspace

In addition to the triangular configuration there is one with all electrons in a plane. For fixed external field the potential (4) possesses a stationary point with one electron on the polarization axis and the other two placed on either side symmetrically with respect to the field axis (Fig. 3). Without loss of generality we may assume that the electrons are confined to the $`x`$-$`z`$-plane and then the $`C_{2v}`$ subspace is spanned by

$`y_1`$ $`=`$ $`y_2=y_3=0`$ (22)
$`p_{y_1}`$ $`=`$ $`p_{y_2}=p_{y_3}=0,`$ (23)
$`x_1`$ $`=`$ $`0,p_{x_1}=0,`$ (24)
$`x_2`$ $`=`$ $`x_3=x,p_{x_2}=p_{x_3}=p_x/2,`$ (25)
$`z_2`$ $`=`$ $`z_3=z,p_{z_2}=p_{z_3}=p_z/2.`$ (26)
The phase space for motion in this subspace is six-dimensional, with variables $`(x,z,z_1,p_x,p_z,p_{z1})`$. The Hamiltonian becomes

$`H`$ $`=`$ $`{\displaystyle \frac{p_x^2+p_z^2}{4}}+{\displaystyle \frac{p_{z_1}^2}{2}}{\displaystyle \frac{6}{\sqrt{x^2+z^2}}}{\displaystyle \frac{3}{z_1}}`$ (29)
$`+{\displaystyle \frac{1}{2|x|}}+{\displaystyle \frac{2}{\sqrt{x^2+(zz_1)^2}}}`$
$`(2z+z_1)f(t).`$
In contrast to the $`C_{3v}`$ subspace the saddle cannot be given analytically and has to be found numerically. For $`f(t)=1`$, the saddle is located at $`|x_s|=1.1607`$, $`z_s=1.1143`$ and $`z_{1,s}=1.4665`$ and has the potential energy $`V_s=7.3902`$. The configuration is elongated along a line perpendicular to the field axis, with axial position of the central electron further out than that of the outer two. Nevertheless, the distance of the central electron to the nucleus is smaller than that of the outer electrons. It thus experiences a stronger attractive force and the balance of force requires that the other two electrons contribute a component pointing away from the nucleus, i.e. that $`|z_s|<|z_{1,s}|`$, as indeed observed.

Harmonic approximation in the $`C_{2v}`$ subspace around the saddle reveals one stable and two unstable modes. One of the unstable directions, with Lyapunov exponent $`\lambda _{2r}=1.0980`$, corresponds to simultaneous escape of three electrons: in the eigenspace

$`z_1`$ $`=`$ $`z_{1,s}+w_{2r}`$ (30)
$`x`$ $`=`$ $`x_s+0.6183w_{2r}`$ (31)
$`z`$ $`=`$ $`z_s+1.1417w_{2r},`$ (32)
all three components increase for positive $`w_{2r}`$. This eigenspace thus corresponds to the reaction coordinate in the triangular configuration. The other eigenmode

$`z_1`$ $`=`$ $`z_{1,s}+u_{2a}`$ (33)
$`x`$ $`=`$ $`x_s0.2313u_{2a}`$ (34)
$`z`$ $`=`$ $`z_s0.3127u_{2a},`$ (35)
with Lyapunov exponent $`\nu _{2a}=1.7937`$ is related to double or single ionization, according to the sign of $`u_{2a}`$.

Harmonic approximation in the full space gives in addition to the modes within the subspace three stable modes, two unstable ones reflecting divergence from the $`C_{2v}`$ subspace and one neutral mode connected with an overall rotation around the field axis. The neutral mode is connected with conservation of angular momentum and was eliminated in the $`C_{3v}`$ case by transformation to polar coordinates; it shows up in the $`C_{2v}`$ analysis where cartesian coordinates are more convenient since one electron would come to lie on the singularity of the polar coordinate system. The first unstable mode with Lyapunov exponent $`\nu _{2b}=0.9024`$ corresponds to deviations along the $`y`$-axis,

$`y_1`$ $`=`$ $`u_{2b}`$ (36)
$`y_2`$ $`=`$ $`0.5752u_{2b}`$ (37)
$`y_3`$ $`=`$ $`0.5752u_{2b}.`$ (38)
The other unstable mode, with the Lyapunov exponent $`\nu _{2c}=1.3712`$, is related to symmetry breaking within the $`x`$-$`z`$-plane,

$`\begin{array}{cc}x_1=0.0746u_{2c}& z_1=z_{1,s}\\ x_2=x_s+0.6598u_{2c}& z_2=z_s+u_{2c}\\ x_3=x_s+0.6598u_{2c}& z_3=z_su_{2c}.\end{array}`$ (39)
The Lyapunov exponents are of similar magnitude as in the $`C_{3v}`$ case and the adiabatic approximation remains justified.

### C Comparison between the subspaces

Analysis of the stationary points in phase space shows that there are two pathways to non-sequential triple ionization. The triangular configuration has a potential barrier (for scaled field $`Ff(t)=1`$) of $`7.6673`$, slightly lower than the one $`7.3902`$ for the saddle in the $`C_{2v}`$ subspace. Therefore, for increasing energy in the compound complex, triple ionization via the triangular path becomes possible first. The triangular state is also less unstable than the $`C_{2v}`$ configuration: it has a largest Lyapunov exponent for symmetry breaking perturbations of $`\nu _a=1.4496`$, as compared to $`\nu _{2a}=1.7937`$. For motion along the reaction coordinates the Lyapunov exponent of the triangular configuration is $`\lambda _r=1.1054`$ and thus slightly larger than $`\lambda _{2r}=1.0980`$ for the $`C_{2v}`$ configuration. The relative phase space weight of trajectories near the symmetry subspaces is determined by a competition between motion along the reaction coordinate and amplification of deviations from symmetry . Generally, the faster the motion along the reaction coordinate and the slower the symmetry breaking the larger the phase space that is dominated by the saddle.

There are no quantitative estimates for the phase space region influenced by the saddle, except close to a threshold, where a suitable extension of Wanniers arguments can be employed . For triple ionization and in the presence of several symmetry breaking modes, the generalization  gives a threshold behavior $`\sigma (\mathrm{\Delta }E)(\mathrm{\Delta }E)^\alpha `$ with exponent

$$\alpha =\left(\nu _i\right)/\lambda _j$$
(40)
where $`\lambda _j`$ is the Lyapunov exponent of the reaction coordinate and the $`\nu _i`$ are the ones of the symmetry breaking modes. Specifically, for our example here we find for the triangular configuration

$$\alpha _3=2.6228$$
(41)
and for the $`C_{2v}`$ configuration

$$\alpha _2=3.7043.$$
(42)
This already suggests that the triangular configuration will be the more dominant configuration. Besides the threshold behavior the two saddles also differ considerably in the final momentum distributions of ionizing electrons and this will be analyzed in the next section.

## III Ion momenta distributions

In the present section we show numerical simulations of the non-sequential ionization. We restrict our calculations to the symmetry subspaces and compare the resulting final ion momentum distributions. Even though the symmetry subspaces are of zero phase space weight and are furthermore unstable non-sequential triple ionization trajectories in full phase space have to pass sufficiently close to the saddles and these subspaces so that they will show distributions very close to the ones obtained in the subspaces . We can, therefore, also compare our data with experimental distributions. In order to allow comparison with experiments we now omit the elimination of the maximal field amplitude $`F`$ by the rescaling (10).

In the experiments of Moshammer et al.  on triple ionization of Ne ultrashort (30 fs) laser pulses at 795 nm wavelength and with peak intensities of 1.5 PW/cm2 were used. This corresponds to a frequency of $`\omega =0.057`$ a.u. (atomic units), a pulse duration, measured as full width at half maximum, of 11 periods $`2\pi /\omega `$ and a maximal field strength of $`F=0.207`$ a.u.. Among the measured data we focus on the distributions of ion momenta parallel to the polarization axis. In the limit of negligible small momentum transfer by the absorbed photons, the ion momentum $`\stackrel{}{p}_{ion}`$ reflects the sum of the momenta of the emitted electrons, $`\stackrel{}{p}_{ion}=_i\stackrel{}{p}_i`$ .

The three-electron Hamiltonian (1) corresponds to lithium. In order to relate the calculations to other atoms with more electrons interactions with core electrons have to be neglected and the point of reference in energy has to be shifted to the threshold for triple ionization. Specifically, for the modeling of the experiments  on triple ionization in Ne we thus assume that in a rescattering process the energy transfer is less than the threshold for triple ionization (about $`4.6`$ a.u.). The precise value of the energy transfered in the rescattering event and thus of the energy $`E`$ of the initial compound state cannot be determined within our model and constitutes a free parameter. Similarly, the time $`t_0`$ during the pulse when the rescattering complex forms is a second free parameter. But both parameters can be determined rather reliably by comparison with experiments . We now turn to the simulations within the two symmetry subspaces.

### A Simulations in the C3v subspace

In a first series of simulations we focus on the triangular configuration. We fix the initial energy $`E`$ and analyze the dependence of the results on the initial time $`t_0`$. To this end we chose a microcanonical distribution of initial conditions in the $`C_{3v}`$ symmetry subspace for energy $`E=0.5`$ a.u. and random phases $`\varphi `$ for different initial time $`t_0`$. The final distributions of ion momenta parallel to the polarization axis are shown in Fig. 4. As in our previous analysis of the non-sequential double ionization process  the distributions has a double hump structure and the width of the distributions increases when the initial time $`t_0`$ approaches the peak amplitude of the field. The width of the distributions can be estimated from the maximal energy a free electron can acquire in the field, i.e. twice the ponderomotive energy. For three electrons ejected in the same direction, the corresponding maximal ion momentum depends on the amplitude of the field at the point in time when they are ionized, i.e.

$$p_{max}=\frac{3Ff_p(t)}{\omega }$$
(43)
where $`f_p(t)=\mathrm{sin}^2(\pi t/T_d)`$ is the pulse envelop. If this time is taken to equal the starting time $`t_0`$ of the simulations we find values of $`p_{max}`$ which correspond very well to the widths of the distributions in Fig. 4. Fig. 4 indicates also that the distributions are basically the same independent of whether $`t_0`$ is chosen before or after the maximum of the pulse, provided $`f_p(t_0)`$ is the same. This implies that the dominant ionization takes place during the first field cycle after formation of the complex. Fitting the width of the calculated distribution to the experimental results allows one to estimate the moment in the pulse when the majority of the triply ionized ions are created.

The dependence of the ion momentum distributions on the initial energy $`E`$ for fixed initial time $`t_0`$ is presented in Fig. 5. The widths of the distributions does not change significantly with $`E`$ but their shape does. For initial energy close to the minimal saddle energy, which for $`F=0.207`$ a.u. is $`3.49`$ a.u., the distributions show a single maximum only while for higher energies two humps form. The electrons that cross the saddle when the energy is near the saddle energy are slow and the combined interactions of the external and Coulomb fields shapes the distribution . The situation is different for the high energy case: then shortly after the electrons cross the saddle the interaction with the electric field is stronger than the attraction to the nucleus and the distribution is mostly shaped by the laser field . As the initial energy of the complex is determined in a rescattering event and is higher for stronger external fields, the distribution of ion momenta should show a transition from a distribution with a single maximum near the threshold for triple ionization to one with a double maximum higher up. For even higher fields the constraints from the triangular configuration are relaxed and less symmetric modes of triple ionization become possible.

### B Numerical simulation in the $`C_{2v}`$ subspace

We now turn to the ion momentum distributions in the $`C_{2v}`$ subspace. Much of the analysis proceeds as in the previous subsection, except for the difference in Hamiltonians (here it is (29)) and the choice of initial conditions. The microcanonical ensemble of the initial conditions is not straightforward to realize because (even for fixed time) for nonzero external field the system is open. In the previous section, the difficulties could be overcome by constraining initial conditions to lie on the energy shell and in the hypersurface $`Z=0`$ (for details see ). Here such a restriction is not sufficient, phase space remains open. The troublesome configurations are those where electron distances to the nucleus are strongly asymmetric, one being close and two far away or vice versa. Since we assume that the initial state is formed during a rescattering and should be confined to region close to the core, such configurations cannot be formed. We, therefore, choose initial conditions microcanonically and require in addition that distances of the electrons to the nucleus are not larger than the minimal distance of the saddle at maximal field to the nucleus.

As in the previous section, the width of the momentum distributions determines the initial time $`t_0`$ when triple ionization takes place. The dependence on initial energy $`E`$ is also similar, with a collapse of all structures into a single peak near zero momentum as the initial energy decreases towards the energy of the saddle (Fig. 6).

In Fig. 6 we show final parallel ion momentum distributions corresponding to non-sequential triple ionization in the $`C_{2v}`$ subspace for three different initial energy $`E`$. The structure of the distributions is significantly different from that in the previous section. In Fig. 6(a) and 6(b), besides two local maxima on the edges of distributions, a strong contribution zero momentum value that is split into two maxima appears. The origin of these structures can be related to qualitatively different electron paths. Initially, after crossing the saddle all three electrons escape towards the same side of the nucleus. However, the final momenta along the polarization axis may point to either the same or opposite directions, depending on the initial conditions and the phase of the field, as shown in Fig. 7. If the final momenta point in the same direction, they can be much larger than if they point in opposite directions. Therefore, the edge maxima are the effect of the electrons escape towards the same direction while the central structure corresponds to an escape in opposite directions.

### C Comparison with experiment

Given the qualitatively very different distributions of final ion momenta in the two subspaces we can use a comparison with experimental data  to draw conclusions about the dominant ionization mode. In Fig. 8 we show the experimental parallel ion momentum distribution together with the one calculated from the $`C_{3v}`$ symmetry subspace. In the calculation we have chosen the initial time $`t_0=0.33T_d`$ that gives the width of the resulting distribution in agreement with the experiment. The initial energy has been taken as $`E=1`$ a.u. but it can be taken in a wide range of values above $`1.5`$ a.u. without significant changes in momentum distributions (Fig. 5). The main features of the experimental distribution are well reproduced in our calculations. However, the minimum in the distribution is less pronounced than that observed in the experiment. Overall, we take this good agreement as strong indication that triple ionization occurs in the neighborhood of the symmetric process discussed here.

The distribution for perpendicular ion momenta, which have been measured in the experiment , can not be calculated within the symmetry subspace since this component vanishes exactly by symmetry. In full phase space ionizing electrons pass close to the saddle but not exactly symmetrically and consequently the transverse ion momenta do not vanish.

## IV Conclusion

We have considered the pathways along which a highly excited three electron complex can decay towards triple ionization. In non-sequential multiple ionization of atoms in strong laser fields such a complex forms when the rescattered electron collides with the core. The subsequent non-sequential decay then has to follow the pathways discussed here. We found two possible pathways, one proceeding via a symmetric triangular configuration, and one via a planar elongated configuration. A stability analysis shows that the ionization close to the $`C_{3v}`$ symmetry subspace occurs at lower field intensities and has a larger cross section (close to threshold). Trajectory simulations within this subspace give the final ion momentum distributions in very good agreement with the recent experimental data . The $`C_{2v}`$ configuration has a higher threshold and smaller cross section and gives a momentum distribution different from the experimentally observed one, so that we conclude that the triangular configuration is indeed the dominant one.

We have not discussed pathways to sequential triple ionization, and one can expect that there are several, including sequential ionization of one electron after the other, a double ionization followed by a single ionization or the other way around, or various combinations of partial ionizations and rescatterings. As in the case of double ionization  one can also have sequential ionization from failed attempts to triple ionize along one of the symmetric saddles. The experimental results of  suggest that, at least in the intensity range considered here, correlated non-sequential electron escape dominates. This is further supported by our previous analysis of double ionization , where we found that the double hump structure observed in the experimental ion momentum distributions  does not appear in the sequential process. That is, the electron momenta are not correlated and the resulting distributions show a strong single maximum instead of the double hump structure.

Finally, let us mention that part of our analysis can be easily extended to non-sequential multiple ionization where the number of escaping electrons, $`N`$, is greater than three . The triangular configuration extends to a placement of electrons on a ring with $`C_{Nv}`$ in plane perpendicular to the field axis. Calculation of the positions and the energy of the saddle shows that the energy is non-monotonic and begins to increase with $`N`$ for more than 9 electrons. For more than 13 electrons the saddle does not exist any more. Therefore, for many electrons the configurations with highest symmetry will disappear and others will dominate the ionization process. For triple ionization we have identified one, the planar configuration, and for more than three electrons we expect a rapidly increasing number of stationary configurations. The identification of a dominant pathway will then require an appropriate stability analysis and an application of the generalized Wannier threshold law .

## V Acknowledgment

We would like to thank H. Giessen, W. Becker, R. Kopold and H. Rottke for discussions. Financial support by the Alexander von Humboldt Foundation and KBN under project 5 P03B 088 21 are gratefully acknowledged.


text:
# Electrorotation of a pair of spherical particles

## I Introduction

When a suspension of colloidal particles or biological cells is exposed to an external electric field, the analysis of the frequency-dependent response yields valuable information on various processes, like the structural (Maxwell-Wagner) polarization effects . For most substances, the permittivity and conductivity are only constant over a limited frequency range. The general tendency is for permittivity to decrease, and conductivity to concomitantly increase, in a series of steps as frequency increases. These step changes are called dispersions.

While the frequency-dependent response of biological cells can be investigated by the method of dielectric spectroscopy , conventional dielectrophoresis and electrorotation (ER) analyze the frequency dependence of translations and rotations of single particles in an inhomogeneous and rotating external field, respectively . With the recent advent of experimental techniques such as automated video analysis  as well as light scattering methods , the motion of particles can be accurately monitored. In ER, the rotating field induces a cell dipole moment which rotates at the angular frequency of the external field. Any dispersion process causes a spatial phase shift between the induced dipole moment and the external field vector, giving rising to a torque which causes the rotation of individual cells .

In the dilute limit in which a small volume fraction of particles or cells are suspended in a medium while the particles are well separated, we can concentrate on the ER response of individual particles by ignoring the mutual interactions between the particles. However, if the suspension is non-dilute, one may not ignore the interactions. Moreover, when the strength of the rotating electric field increases, the Brownian motion can be ignored and the polarized particles tend to aggregate along the rotating field even in the dilute limit.

As an initial model, we will consider a pair of interacting particles dispersed in a suspension. When the two particles approach and finally touch, the mutual polarization interaction between the particles leads to a change in the dipole moment of individual particles and hence the ER spectrum, as compared to that of the well-separated particles. We will present a theoretical study on the ER spectrum of two spherical particles in the presence of a rotating electric field, and capture the mutual polarization effects by the method of multiple images , recently developed to calculate the interparticle force between touching spherical particles.

## II Multiple image method for a pair of dielectric spheres

We first consider an isolated spherical cell of complex permittivity $`\stackrel{~}{ϵ}_1=ϵ_1+\sigma _1/(i\omega )`$ dispersed in a suspension medium of $`\stackrel{~}{ϵ}_2=ϵ_2+\sigma _2/(i\omega )`$, where $`\omega `$ is the frequency of the external field $`\stackrel{}{E}_0`$, and $`i=\sqrt{1}`$. If one applies a rotating electric field, there is in general a phase difference between the induced dipole moment and the applied field vector, and thus a net torque is acting on the cell. The spin friction (analogous to the Stokes drag for translational motions) causes a steady-state rotational motion. In this case, the dipole moment of the isolated cell $`\stackrel{~}{p}`$ is

$$\stackrel{~}{p}=\frac{1}{8}\stackrel{~}{ϵ}_2D^3bE_0,$$
(1)
where $`b=(\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_2)/(\stackrel{~}{ϵ}_1+2\stackrel{~}{ϵ}_2)`$ is the dipole factor, and $`D`$ the diameter of the cell.

To account for the effect of multipole interaction on the electrorotation, we consider a pair of touching spherical cells at a separation $`R`$ suspended in a medium. We will calculate the effect of multipole interaction on the dipole moment. For a rotating field, the total dipole moment of one cell of the pair is given by

$$\stackrel{~}{p}^{}=\stackrel{~}{P}_T\mathrm{cos}^2\theta +\stackrel{~}{P}_L\mathrm{sin}^2\theta =\frac{1}{2}(\stackrel{~}{P}_T+\stackrel{~}{P}_L),$$
(2)
where $`\theta `$ is the angle between the dipole moment and the line joining the centers of the spheres. The average over all possible orientations of the cells gives a factor 1/2. In Eq.(2), $`\stackrel{~}{P}_T`$ ($`\stackrel{~}{P}_L`$) is the transverse (longitudinal) dipole moment, being perpendicular (parallel) to the line joining the centers of the spheres. The transverse and longitudinal dipole moments are given by

$`\stackrel{~}{P}_L`$ $`=`$ $`\stackrel{~}{p}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(2b)^n\left({\displaystyle \frac{\mathrm{sinh}\alpha }{\mathrm{sinh}(n+1)\alpha }}\right)^3,`$ (3)
$`\stackrel{~}{P}_T`$ $`=`$ $`\stackrel{~}{p}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(b)^n\left({\displaystyle \frac{\mathrm{sinh}\alpha }{\mathrm{sinh}(n+1)\alpha }}\right)^3,`$ (4)
where $`\alpha `$ satisfies the relation $`\mathrm{cosh}\alpha =R/D`$. For two touching particles, the dipole moment of a particle is given by

$$\stackrel{~}{p}^{}=\frac{1}{8}\stackrel{~}{ϵ}_2D^3b^{}E_0,$$
(5)
where the dipole factor of a pair is given by:

$$b^{}=\frac{8\stackrel{~}{p}^{}}{\stackrel{~}{ϵ}_2D^3E_0}=\frac{1}{2}b\left[\underset{n=0}{\overset{\mathrm{}}{}}(2b)^n\left(\frac{\mathrm{sinh}\alpha }{\mathrm{sinh}(n+1)\alpha }\right)^3+\underset{n=0}{\overset{\mathrm{}}{}}(b)^n\left(\frac{\mathrm{sinh}\alpha }{\mathrm{sinh}(n+1)\alpha }\right)^3\right].$$
(6)
It is known that the angular velocity of electrorotation of a particle is related to the dipole factor ($`b`$ or $`b^{}`$) as follows

$$\mathrm{\Omega }=F(ϵ_2,\eta ,E_0)Im[DipoleFactor]$$
(7)
where $`F`$ is a function of $`ϵ_2`$, the dynamic viscosity $`\eta `$ of the medium, as well as the electric-field magnitude $`E_0`$, and $`Im[\mathrm{}]`$ denotes the imaginary part of $`[\mathrm{}]`$. For a spherical particle, $`F(ϵ_2,\eta ,E_0)=ϵ_2E_0^2/2\eta `$. Regarding $`\eta `$, the spin friction expression suffices for an isolated particle. However, for two interacting particles, we must consider the more complicated suspension hydrodynamics.

## III Spectral representation and dispersion spectrum

In a recent paper , we studied the dielectric behavior of cell suspensions by employing the Bergman-Milton spectral representation of the effective dielectric constant . By means of the spectral representation, we derived the dielectric dispersion spectrum in terms of the electrical and structure parameters of the cell models. The spectral representation is a rigorous mathematical formalism of the effective dielectric constant of a composite material. The essence of the spectral representation is to define the following transformations. If we denote a complex material parameter

$$\stackrel{~}{s}=\left(1\frac{\stackrel{~}{ϵ}_1}{\stackrel{~}{ϵ}_2}\right)^1,$$
(8)
then the dipole factor $`b^{}`$ admits the general form

$$b^{}=\underset{n}{}\frac{F_n}{\stackrel{~}{s}s_n},$$
(9)
where $`n`$ is a positive integer, i.e., $`n=1,2,\mathrm{}`$, and $`F_n`$ and $`s_n`$ are the $`n`$-th microstructure parameters of the composite material .

Thus the spectral representation offers the advantage of the separation of material parameters (namely the dielectric constant and conductivity) from the cell structure information, thus simplifying the study. From the spectral representation, one can readily derive the dielectric dispersion spectrum, with the dispersion strength as well as the characteristic frequency being explicitly expressed in terms of the structure parameters and the materials parameters of the cell suspension . The actual shape of the real and imaginary parts of the permittivity over the relaxation region can be uniquely determined by the Debye relaxation spectrum, parametrized by the characteristic frequencies and the dispersion strengths. So, we can study the impact of these parameters on the dispersion spectrum directly. In what follows, we further express the dipole factor $`b`$ and $`b^{}`$ in the spectral representation. The dipole factor $`b`$ admits the form $`b=F_1/(\stackrel{~}{s}s_1)`$ in the spectral representation, where $`F_1=1/3`$ and $`s_1=1/3`$. For the dielectric dispersion spectrum, we further define two dimensionless real parameters

$$s=(1ϵ_1/ϵ_2)^1,\mathrm{and}t=(1\sigma _1/\sigma _2)^1,$$

where $`s`$ and $`t`$ are the dielectric and conductivity contrasts respectively. After simple manipulations, Eq.(7) becomes

$$b=\frac{F_1}{ss_1}+\frac{\delta ϵ_1}{1+if/f_c}$$
(10)
where $`f_c`$ is the characteristic frequency, at which the maximum angular velocity of electrorotation occurs, and $`\delta ϵ_1`$ is the dispersion magnitude :

$`\delta ϵ_1`$ $`=`$ $`F_1{\displaystyle \frac{st}{(ts_1)(ss_1)}},`$ (11)
$`f_c`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{\sigma _2}{ϵ_2}}{\displaystyle \frac{s(ts_1)}{t(ss_1)}}.`$ (12)
Similarly, the dipole factor of a pair of interacting particles $`b^{}`$ can be rewritten in the spectral representation as

$$b^{}=\underset{m=1}{\overset{\mathrm{}}{}}\left(\frac{F_m^{(T)}}{\stackrel{~}{s}s_m^{(T)}}+\frac{F_m^{(L)}}{\stackrel{~}{s}s_m^{(L)}}\right)$$
(13)
where

$`F_m^{(T)}`$ $`=`$ $`F_m^{(L)}={\displaystyle \frac{2}{3}}m(m+1)\mathrm{sinh}^3\alpha e^{(2m+1)\alpha },`$ (14)
$`s_m^{(T)}`$ $`=`$ $`{\displaystyle \frac{1}{3}}[1+e^{(1+2m)\alpha }],s_m^{(L)}={\displaystyle \frac{1}{3}}[12e^{(1+2m)\alpha }].`$ (15)
In the above derivation, we have used the identity

$$\frac{1}{\mathrm{sinh}^3x}=\underset{m=1}{\overset{\mathrm{}}{}}4m(m+1)e^{(1+2m)x}.$$

We should remark that the spectral representation of $`b^{}`$ \[Eq.(10)\] is an exact transformation  of the multiple image expression \[Eq.(4)\]. Thus, the spectral representation of $`b^{}`$ consists of a discrete set of simple poles (see Fig.5 below). The $`m=1`$ pole $`s_1^{(L)}`$ ($`s_1^{(T)}`$) for the longitudinal (transverse) field case deviates significantly from $`s=1/3`$ while $`s_1^{(L)}`$ becomes small in the touching limit, $`\alpha 0`$. In what follows, we will show that this pole gives a significant contribution to the ER spectrum at low frequency. As $`m`$ increases, however, the series $`s_m^{(T)}`$ and $`s_m^{(L)}`$ accummulate to $`s=1/3`$, giving rise to a dominant contribution near the isolated-sphere frequency $`f_c`$. Moreover, each term in the spectral representation expression of $`b^{}`$ can be rewritten as,

$`{\displaystyle \frac{F_m^{(T)}}{\stackrel{~}{s}s_m^{(T)}}}`$ $`=`$ $`{\displaystyle \frac{F_m^{(T)}}{ss_m^{(T)}}}+{\displaystyle \frac{\delta ϵ_m^{(T)}}{1+if/f_{mc}^{(T)}}},`$ (16)
$`{\displaystyle \frac{F_m^{(L)}}{\stackrel{~}{s}s_m^{(L)}}}`$ $`=`$ $`{\displaystyle \frac{F_m^{(L)}}{ss_m^{(L)}}}+{\displaystyle \frac{\delta ϵ_m^{(L)}}{1+if/f_{mc}^{(L)}}},`$ (17)
where $`\delta ϵ_m^{(T)}`$ (or $`\delta ϵ_m^{(L)}`$) and $`f_{mc}^{(T)}`$ (or $`f_{mc}^{(L)}`$) are the dispersion magnitudes and the characteristic frequencies of the transverse (longitudinal) field cases, obtained respectively by replacing $`F_1`$ and $`s_1`$ in the expressions of $`\delta ϵ_1`$ and $`f_c`$ in Eq.(9) with $`F_m^{(T)}`$ (or $`F_m^{(L)}`$) and $`s_m^{(T)}`$ (or $`s_m^{(L)}`$). Note that $`f_{1c}^{(L)}`$ is significantly lower than $`f_c`$. The ER spectrum of two spheres will consist of a series of sub-dispersions. In what follows, we will investigate the effect of multipole interaction on the electrorotation of cell, as shown in the next section.

## IV Numerical results

For convenience, we let $`F(ϵ_2,\eta ,E_0)=1`$ in our numerical calculation, both for the isolated-particle and touching-particle cases, in order to study the effect of multipole interaction on the ER spectrum. It is because we will focus on the frequency dependence of the ER spectrum. On the other hand, we do not know the dynamic viscosity for two approaching spheres.

Fig.1 is plotted to compare the two cases for three different conductivity contrast $`t`$. It is evident that, for both cases, a small $`|t|`$ yields a high characteristic frequency, at which the peak locates. Moreover, a smaller $`|t|`$ gives a larger peak value. Given a constant $`t`$, the touching-cell case predicts a shift of the peak location to lower frequency than that predicted by the isolated-cell case. At the same time, the peak for the touching-cell case broaden. This behavior arises from the effect of the multipole interaction. In addition, if $`t`$ is given, the two cases also predict different behavior within the low-frequency region, but almost the same behavior within the high-frequency region. Especially for small $`|t|`$ (e.g., $`t=1/900`$), a second peak occurs at lower frequency. This is in accord with the prediction of the spectral representation.

Fig.2 is plotted to exhibit the effect of the material parameter $`s`$ on the two cases. It is observed that, for both cases, a larger $`s`$ gives a larger peak value. Moreover, the peak of different $`s`$ appears at almost the same frequency for both cases. In summary, $`s`$ may affect the peak value, but not the peak location.

Fig.3 is plotted for three different values of the medium conductivity $`\sigma _2`$ at $`s=1.1`$ and $`t=1/90`$. For a larger $`\sigma _2`$, we get a larger characteristic frequency at which the peak occurs. It is evident that $`\sigma _2`$ may affect the peak location, but not the peak value.

In Fig.4, as $`R/D1`$ (e.g., $`R/D=1.0333`$), the deviation between the two cases is evident. In other words, the multipole interaction does play an important role in the ER spectrum and the effect cannot be neglected for touching particles. However, as the separation increases, say, $`R/D=2`$, both cases predicts the same ER spectrum. From the results, we would say that the effect of multipole interaction may be neglected at $`R/D>2`$.

In Fig.5, we plot the spectral representation for the two sets of poles $`s_m^{(L)}`$ (or $`s_m^{(T)}`$) for $`m=1`$ to $`100`$. As $`m`$ increases, $`s_m^{(L)}`$ increases up to $`1/3`$, while $`s_m^{(T)}`$ decreases towards $`1/3`$. For small $`R/D`$ ratio (say, $`R/D=1.0333`$), the longitudinal field plays an more important role in determining electrorotation spectrum than the transverse field. However, if $`R/D`$ ratio is large enough, say, $`R/D2`$, (results not shown here), the effect of multipole interaction is negligible both for the longitudinal and transverse field cases and the results for an isolated cell recover. In summary, as the two cells approach and finally touch, one must take into account the effect of multipole interaction on the electrorotation of cells.

## Discussion and Conclusion

Here a few comments on our results are in order. In this work, we attempted a theoretical study of the electrorotation of two approaching spherical particles in the presence of a rotating electric field. When the two particles approach and finally touch, the mutual polarization interaction between the particles leads to a change in the dipole moment of individual particles and hence the ER spectrum, as compared to that of the well-separated particles. The mutual polarization effects are studied via the multiple image method. From the results, we find that the mutual polarization effects can change the characteristic frequency substantially. We should remarked that the multiple image method for two dielectric spheres is only approximate , but the approximation is quite good . More accurate calculations based on bispherical coordinates can be attempted. However, we believe that similar conclusions can be obtained.

Moreover, when the volume fraction of the suspension is large, the rotating electric field leads to plate-like aggregation of particles in the plane of the applied field. In this case, the electrorotation spectrum can be modified further due to the local field effects arising from the many-particle system. To this end, a first-principles approach can be used to handle the many-particle and multipole interactions .

Regarding the dynamic viscosity, the spin friction expression for an isolated spherical cell suffices. However, for many cells interacting in a suspension, we must consider the hydrodynamics. Regarding the dynamic viscosity for two touching spheres, we must consider the suspension hydrodynamics of two approaching spheres. This is a formidable task and should be a topic for future work.

## Acknowledgments

This work was supported by the Research Grants Council of the Hong Kong SAR Government under grant CUHK 4245/01P. G.Q.G. acknowledges the financial support from a Key Project of the National Natural Science Foundation of China under grant 19834070. K. W. Y. acknowledges useful correspondence with Professor Mike Thorpe regarding the spectral representation of two approaching spheres, and extensive discussion with Professor Tuck Choy on the multiple image method.

Fig.1

Fig.2

Fig.3

Fig.4

Fig.5