text:
# On the Origin of Early Superhumps in WZ Sge-Type Stars

## 1 Introduction

Dwarf novae are a class of cataclysmic variables (CVs), which are close binary systems consisting of a white dwarf and a red dwarf secondary transferring matter via the Roche-lobe overflow. WZ Sge-type dwarf novae (cf. Bailey (1979); Downes, Margon (1981); Patterson et al. (1981); Downes (1990); O’Donoghue et al. (1991); Kato et al. (2001b)) are a small subgroup of dwarf novae characterized by the long ($``$ 10 yr) outburst recurrence time and the large ($``$ 8 mag) outburst amplitude. They are a subclass of SU UMa-type dwarf novae (cf. Warner (1995b)).

The most remarkable signature of WZ Sge-type outbursts is the presence of “early superhumps”111 This feature is also referred to as orbital superhumps (Kato et al., 1996) or outburst orbital hump (Patterson et al., 1998).  during their earliest stage of superoutbursts. Early superhumps have a period extremely close to that of the binary period222 The best-established case is the 2001 superoutburst of WZ Sge (Ishioka et al. (2002) and references therein). In a few other systems (e.g. AL Com: Kato et al. (1996), Patterson et al. (1996), Nogami et al. (1997); EG Cnc: Kato et al. (1997); Matsumoto et al. (1998); Patterson et al. (1998)), the periods of early superhumps have been found to be in good agreement with their quiescent photometric periods (most likely representing orbital periods). , and commonly show double-humped profile (figure 1; see also Kato et al. (1996) and Kato et al. (1998) for detailed discussions), in contrast to ordinary superhumps of SU UMa-type dwarf novae \[see Warner (1995b) for basic properties of SU UMa-type dwarf novae\]. Early superhumps are the most discriminative feature of WZ Sge-type outbursts, and have not been detected in other dwarf novae (see Kato et al. (2001a) and Kato et al. (2001b)). Although several models have been historically proposed \[e.g. enhanced hot spot (Patterson et al., 1981), immature form of superhumps (Kato et al., 1996), jet or thickened disk (Nogami et al., 1997)\], none of them has successfully explained all the observed features of early superhumps.

## 2 Recent Observations of Early Superhumps

During the early stage of the 2001 outburst of WZ Sge (e.g. Ishioka et al. (2002)), several authors (Steeghs et al. (2001); Baba et al. (2001)) detected two-armed arch-like structures on Doppler tomograms(Marsh, Horne, 1988) of HeII and CIII/NIII emission lines. The feature resembled those of ‘spiral structure’ seen on Doppler tomograms of the dwarf nova IP Peg (Steeghs et al. (1997); Steeghs et al. (1998)) in outburst. The spiral structure in Doppler tomograms of IP Peg has been widely believed to represent tidally induced spiral shocks or spiral waves (Steeghs et al., 1997). However, it is not clear whether the same interpretation can apply to a completely different (very low mass-ratio, $`q=M_2/M_1<0.1`$ in contrast to $`q0.5`$ in IP Peg, and very low-mass transfer rate, several orders of magnitudes lower than that of IP Peg) binary parameters of WZ Sge.

At the epoch of the detection of two-armed arch-like structures in WZ Sge, the star showed double-wave early superhumps (Ishioka et al., 2002). Although there have been suggestions that early superhumps in photometry and two-armed structures on Doppler tomograms have the same origin (Ishioka et al. (2002); Baba et al. (2002); Kuulkers et al. (2002)), no promising idea has been proposed to explain these features in the same scheme.

## 3 Tidally Distorted Accretion Disks in Close Binaries

Most recently, an alternative idea has been proposed to explain arch-like structures in the outbursting IP Peg disk, by considering the irradiation of the elevated disk formed by the horizontal convergence ($`div(\text{v})`$) of three-body orbits (Smak, 2001). This interpretation successfully reproduced the arch-like structures, and the observed strong intensity in the region of ($`v_x>0,v_y>0`$), which is difficult to explain within the scheme of conventional tidally induced spiral shocks. Ogilvie (2001) studied tidal distortion of fluid disks, and succeeded in reproducing the basically same features. Ogilvie (2001) showed, in realistic fluid disks, that $`m=2`$ inner vertical resonance play a more important role even in high $`q`$ binaries.

## 4 New Interpretation of Early Superhumps

In low $`q`$ systems like WZ Sge, tidal force of the secondary is so weak on a compact disk that a condition is usually achieved that angular momentum is more effectively transported by viscosity than tidal force of the secondary (Lin, Papaloizou, 1979). However, during the early stage of vigorous outbursts of WZ Sge-type stars, the disk can sufficiently expand (Osaki, 1995). In usual SU UMa-type dwarf novae, the growth of 3:1 resonance, which is responsible of superhumps (Whitehurst (1988); Hirose, Osaki (1990); Lubow (1991)), is rapid enough to effectively truncate the disk at the resonance radius. In low $`q`$ WZ Sge-type stars, the slow growth \[the growth rate being proportional to $`q^2`$ (Lubow, 1991)\] can allow the disk to expand beyond the 3:1 resonance. 333 Some observational implications of WZ Sge-type stars have been proposed to be a consequence of the disk beyond the 3:1 resonance (Kato et al., 1998).  In such a condition, tidal distortion of the outer edge of the accretion disk, as proposed by Smak (2001) and Ogilvie (2001), can effectively work, as shown in figure 3, $`q`$=0.1 case of Ogilvie (2001). This effect is predominant in the outer edge of the accretion disk, which can explain stronger and distorted emission lines in low velocity regions (e.g. Baba et al. (2002)). Since high-excitation emission lines like HeII require a high temperature ($``$ 25000 K), the concentration of emission lines and arch-like patterns in low velocity regions is otherwise difficult to explain (Smak, 2001). The arch-like pattern in Doppler tomograms in WZ Sge is also stronger in the region of ($`v_x>0,v_y>0`$), which also can be naturally explained as a result of stronger tides near the secondary, as in Smak (2001) and Ogilvie (2001).

The maximum phases of early superhumps have been also problematic. Early superhumps in WZ Sge-type dwarf novae consist of two maxima of different amplitudes (figure 1). In WZ Sge, the brighter maximum is around binary phase 0.6–0.7 (Ishioka et al., 2002). While models involving an enhanced hot spot as in Patterson et al. (1981) would require such a maximum corresponding to the orbital phase of a stronger emission feature, the available observation shows a complete phase reversal (figure 2). Such a reversal can be naturally understood if one considers the source of emission lines and stronger continuum is an irradiated surface (Smak (2001); Ogilvie (2001)) of the accretion disk. \[Since CVs are known to emit most of their energy in the UV wavelengths (e.g. Warner (1995a)), re-emission of a small fraction of intercepted UV photons in the optical wavelengths is expected to explain the observed amplitudes of continuum variation (early superhumps)\]. The decrease in amplitude of early superhumps in accordance with the fading of the outburst (e.g. Ishioka et al. (2002)) can be also naturally understood as an effect of shrinkage of the disk, resulting in decreasing tidal effects.

Ishioka et al. (2002) also showed that the exact period of early superhumps is very slightly (0.05%) shorter than the orbital period. This effect was observed as a phase-shift of maxima of early superhumps. This change looks likely to correspond to a counterclockwise rotation of the patterns on Doppler tomograms at different epochs (figure 2). The rotation of patterns may reflect the variation of the direction of the ‘maximum height distortion’ toward a smaller disk radius (see figure 2 of Ogilvie (2001)). Further detailed calculations of a time-evolution of the outbursting WZ Sge disk would provide a stringest test to this interpretation.

Most recently, Osaki, Meyer (2001) explained early superhumps (which they call early humps) by considering the manifestation of the tidal 2:1 resonance in accretion disks of binary systems with extremely low mass ratios. The present interpretation is different from that of Osaki, Meyer (2001) in its basic mechanism: Osaki, Meyer (2001) requires tidal dissipation while the present interpretation only requires tidal distortion of a steady flow (hence it does not require a shock). The present interpretation has an advantage over all existing models in its simplicity of the assumed physics and in that it can also explain two-armed distortion patterns on Doppler tomograms as well as photometric signals.

The reason why such features are only observed in WZ Sge-type stars can be explained as a combination of the large stored mass in quiescence (Osaki, 1995) and the delayed exertion of the 3:1 resonance, which is only reasonably achieved in extremely low $`q`$ and low mass-transfer rate (and presumably low quiescent viscosity) systems as in WZ Sge-type stars. Since three-dimensional fluid calculation by Ogilvie (2001) showed that $`m=2`$ inner vertical resonance can effectively work in relatively high $`q`$ systems, we can expect similar phenomena expected in Lin, Papaloizou (1979) in larger $`q`$ systems, which cannot hold particle orbits of the 2:1 resonance within the Roche lobe. Although Osaki, Meyer (2001) also presented an idea of distinguishing WZ Sge-type stars from other CVs based on resonances involved, the present interpretation is different from that of Osaki, Meyer (2001) in that Osaki, Meyer (2001) only considered resonances in Paczynski type orbits (e.g. Ogilvie (2001)), while the present interpretation involves the three-dimensional effect as presented by Ogilvie (2001). The appearance of early superhumps and the reproduction of WZ Sge-type outburst in RZ Leo (Ishioka et al., 2001), a binary with $`q=0.14`$, which is approximately twice larger than those of other known WZ Sge-type stars, can be understood as one of the most distinct consequences between the present interpretation and that of Osaki, Meyer (2001).

The author is grateful to H. Baba, R. Ishioka and D. Steeghs for allowing the use of their figures. The author is also grateful to Y. Osaki for drawing attention to the paper by Lin, Papaloizou (1979). This work is partly supported by a grant-in aid (13640239) from the Japanese Ministry of Education, Culture, Sports, Science and Technology.


format: text

text:
# Formation of gas disks in merging galaxies

## 1 INTRODUCTION

Gas falling into the nuclei of interacting and merging galaxies has attracted much attention. Nuclear inflows fuel the bursts of star formation responsible for the anomalous colors (Toomre & Toomre 1972; Larson & Tinsley 1978) and extraordinary IR luminosities (Joseph & Wright 1985; Sanders & Mirabel 1996) of merging galaxies; the stars formed in this manner may play a crucial role in transforming merger remnants into elliptical galaxies (Kormendy & Sanders 1992). Numerical simulations show that these inflows occur when gravitational torques remove angular momentum from shocked galactic gas (Noguchi 1988; Hernquist 1989; Combes, Dupraz, & Gerin 1990; Barnes & Hernquist 1991, 1996). In typical cases $`60`$ percent of the gas initially distributed throughout the progenitor disks may wind up in a nuclear cloud with dimensions of $`0.1\mathrm{kpc}`$.

Less attention has been paid to the gas which fails to reach the nuclei of merger remnants. Extended HI, HII, and X-ray emission is often seen in merging galaxies, but the gas responsible is usually too dilute to directly fuel violent star formation or AGN activity. There has been little numerical work on the fate of extended gas in mergers, in part because it’s expensive to continue a simulation once a compact nuclear gas cloud has formed.

There are reasons to think that extra-nuclear gas may play a role in the evolution of merger remnants. The well-studied merger remnant NGC 7252 exhibits extended emission from ionized gas with remarkably complex kinematics (Schweizer 1982, 1998): an inner disk with radius $`3\mathrm{kpc}`$ and a well-defined axis of rotation is surrounded by a region which shows rapid minor-axis rotation and an apparent velocity reversal. This inner disk, seen as a vivid spiral in HST images (Whitmore et al. 1993), contains $`10^{9.5}\mathrm{M}_{}`$ of molecular and ionized gas (Dupraz et al. 1990). A comparable amount of atomic gas lingers in the tidal tails; these return gas to the main body of the remnant at at rate of $`4\mathrm{M}_{}\mathrm{yr}^1`$ (Hibbard et al. 1994; Hibbard & Mihos 1995). Whatever the fate of this gas, it’s worth asking how a system with such complex kinematics could form in the first place.

In this paper I use simulations to explore the possibility that extended gas disks and associated structures can form as part of the merger process. Each simulation begins with a pair of disk galaxies falling together on a parabolic relative orbit, and follows their collision and the eventual formation and initial relaxation of a merger remnant. The interstellar material is modeled by including a component obeying the equations of motion of an isothermal gas. This approach is deliberately simple-minded when compared to simulations which include star formation, feedback, and other processes (e.g. Katz 1992; Mihos & Hernquist 1994; Gerritsen & Icke 1997; Springel 2000). Star formation is clearly important for the evolution and eventual fate of merging galaxies, but there’s no evidence that stellar processes are essential for the formation of gas disks like the one in NGC 7252. My intention is to present a simple model which appears capable of explaining the disks observed in merger remnants.

## 2 SIMULATIONS

Isothermal simulations are relatively cheap, so an encounter survey can be performed with modest computational resources. For ease of comparison, I adopted the same set of four encounter geometries used in previous studies (e.g. Barnes 1992; Barnes & Hernquist 1996). Table 1 lists disk inclinations $`i`$ and pericentric arguments $`\omega `$ for the passages used in this paper. Each combination was used in three different encounters: first, a close passage of galaxies with a 1:1 mass ratio; second, a more distant passage, again with a 1:1 ratio; and third, a close passage with a 3:1 mass ratio. Thus there are $`12`$ simulations in this survey. For convenience, I designate each simulation by specifying the disk geometry (from Table 1), the mass ratio (either 1:1 or 3:1), and optionally the pericentric separation (C or D). Thus DIR 1:1 C specifies a direct, close passage of equal-mass galaxies, DIR 1:1 D is a more distant version of the same encounter, and INC 3:1 is an inclined passage of galaxies with a 3:1 mass ratio.

The galaxy models in these experiments are the same as those used in earlier studies (Barnes 1998; Bendo & Barnes 2000). Briefly, each model contains a bulge with a shallow cusp (Hernquist 1990), an exponential disk with constant scale height (Freeman 1970; Spitzer 1942), and a dark halo with a constant-density core (Dehnen 1993; Tremaine et al. 1994). Density profiles for these components are

$$\begin{array}{c}\rho _{\mathrm{bulge}}r^1(r+a_{\mathrm{bulge}})^3,\hfill \\ \rho _{\mathrm{disk}}\mathrm{exp}(R/R_{\mathrm{disk}})\mathrm{sech}^2(z/z_{\mathrm{disk}}),\hfill \\ \rho _{\mathrm{halo}}(r+a_{\mathrm{halo}})^4,\hfill \end{array}$$
(1)
where $`r`$ is the spherical radius, $`R`$ is the cylindrical radius in the disk plane, and $`z`$ is the distance from the disk plane. The gas is initially distributed like the disk component and amounts to $`12.5`$ percent of the disk mass. Initial velocities and velocity dispersions for the collisionless components were determined from the Jeans equations (e.g. Hernquist 1993); the disk gas was set rotating at the local circular velocity.

Unless otherwise indicated, the results below are given in simulation units with $`G=1`$. In these units, the galaxies in the equal-mass encounters had total masses $`M_1=M_2=M_{\mathrm{bulge}}+M_{\mathrm{disk}}+M_{\mathrm{halo}}=\frac{1}{16}+\frac{3}{16}+1=\frac{5}{4}`$ and length scales $`a_{\mathrm{bulge}}=0.04168`$, $`R_{\mathrm{disk}}=0.08333`$, $`z_{\mathrm{disk}}=0.007`$, and $`a_{\mathrm{halo}}=0.1`$. These choices yield a model with half-mass radius $`r_{\mathrm{half}}0.28`$, rotation period $`t_{\mathrm{rot}}(r_{\mathrm{half}})1.2`$, and binding energy $`E1.07`$. The same model was used for the larger galaxy in the simulations with a 3:1 mass ratio, while the smaller model was simply scaled down by a factor of $`3`$ in mass and a factor of $`\sqrt{3}`$ in radius, thereby following a $`Mv^4`$ relation similar to the observed luminosity-rotation velocity relation for disk galaxies. The large galaxy model may be roughly scaled to the Milky Way by equating the simulation units of length, mass, and time to $`40\mathrm{kpc}`$, $`2.2\times 10^{11}\mathrm{M}_{}`$, and $`2.5\times 10^8\mathrm{yr}`$, respectively. Gravitational forces were calculated after smoothing the mass distribution using a Plummer kernel with $`ϵ=0.0125`$.

The gas had an isothermal equation of state, $`P=c_\mathrm{s}^2\rho `$, where the sound speed is fixed at $`c_\mathrm{s}=0.0966`$ velocity units. The simulation code evaluates the heating $`\dot{u}`$ of each gas particle due to adiabatic compression and shocks; this energy is assumed to be radiated instantly, so the gas stays at a constant temperature. The sound speed $`c_\mathrm{s}`$ is an order of magnitude smaller than typical circular velocities in these models; thus gas pressure forces are relatively small, and the gas travels on roughly ballistic trajectories except where diverted by shocks. If the large galaxy model is scaled to the MW, the sound speed is $`c_\mathrm{s}15.5\mathrm{km}\mathrm{s}^1`$, corresponding to a gas temperature of $`20000\mathrm{K}`$. I chose this value to take account of non-thermal pressures due to magnetic fields and turbulence, but the more conventional choice of $`10^4\mathrm{K}`$ would yield nearly identical results.

Initial conditions for the encounters were generated by building pairs of galaxy models, placing them on the chosen orbits, and rotating them to the desired orientations. The initial center-of-mass position and velocity of each galaxy were determined from a parabolic orbit of two point masses $`M_1`$ and $`M_2`$ with pericentric separation $`r_{\mathrm{peri}}`$ and time of pericenter $`t_{\mathrm{peri}}`$. For the 1:1 encounters the close passages had $`r_{\mathrm{peri}}=0.2`$, while the distant passage had $`r_{\mathrm{peri}}=0.4`$. The 3:1 encounters used $`r_{\mathrm{peri}}=0.2`$. In all cases $`t_{\mathrm{peri}}=1`$, so point masses on these orbits would reach pericenter exactly one time unit after the start of the simulation. Of course, the actual trajectories deviate from the Keplerian ideal as soon as the galaxy models begin to interpenetrate (e.g. Barnes 1988, 1992). At first passage, these deviations are modest – the close encounters come within $`r_{\mathrm{peri}}0.25`$ at $`t_{\mathrm{peri}}1.03`$, while the distant ones reach $`r_{\mathrm{peri}}0.45`$ at $`t_{\mathrm{peri}}1.04`$.

Each simulation used $`N_{\mathrm{gas}}+N_{\mathrm{star}}+N_{\mathrm{halo}}=24576+29696+32768=87040`$ particles. Collisionless components were followed using standard N-body techniques, while the gas was simulated using Smoothed Particle Hydrodynamics or SPH (e.g. Monaghan 1992); the code uses a hierarchical algorithm to compute gravitational forces, adaptive smoothing to resolve a fixed mass scale in the gas, and adaptive time-steps determined by a Courant condition. To take maximum advantage of a cluster of eight $`333\mathrm{MHz}`$ processors, I ran multiple calculations simultaneously. Most runs required between $`400`$ and $`800`$ processor hours, largely spent in the later stages where high gas densities demand very short time-steps. Energy and angular momentum were conserved to a fraction of a percent in all cases.

The resolution of these simulations is comparable to that available in earlier SPH studies of dissipative encounters involving disk galaxies (e.g. Hernquist 1989; Barnes & Hernquist 1991, 1996; Mihos & Hernquist 1994, 1996). Gravitational forces were evaluated with a fixed smoothing of $`ϵ=0.0125`$ length units, which is small enough to follow the overall dynamics of the interacting galaxies; the nuclei of the resulting merger remnants are poorly resolved, but the disks formed in these simulations are an order of magnitude larger than $`ϵ`$ and thus not severely compromised by gravitational smoothing. Hydrodynamic forces were calculated by smoothing over $`40`$ gas particles; each gas particle has a mass of $`2^{19}1.9\times 10^6`$ mass units, so the SPH calculation has a mass resolution of $`7.6\times 10^5`$. For gas with a density $`\rho `$ and sound speed $`c_\mathrm{s}=0.0966`$, the Jeans length is $`\lambda _\mathrm{J}0.1712\rho ^{1/2}`$ and the Jeans mass is $`M_\mathrm{J}0.00263\rho ^{1/2}`$; the gravitational force calculation resolves $`\lambda _\mathrm{J}`$ if $`\rho 190`$, while the SPH calculation resolves $`M_\mathrm{J}`$ if $`\rho 1200`$. Thus the collapse of small-scale structure in the gas is limited by the spatial resolution of the gravitational force calculation, with hydrodynamic resolution playing only a secondary role. Suppression of small-scale collapse, which would be fatal for studies of fragmentation via Jeans instabilities (e.g. Bate & Burkert 1997), may even benefit these simulations by preserving the smooth structure of the gas. The code does not suppress the Jeans instability on scales $`\lambda _\mathrm{J}>ϵ`$; on such scales, the gas is usually stabilized by the combined effects of pressure and rotation (Toomre 1964).

Nonetheless, these simulations are in fact quite crude. The SPH smoothing length often exceeds the vertical scale height of the simulated gas disks; in effect, these structures are resolved in the radial direction, but not in the vertical direction. Under such circumstances the calculation provides only a rough approximation to the dynamics of a smooth gas, and it seems better to view the SPH code as a locally momentum-conserving scheme in which particles representing gas tend to seek out closed, non-intersecting orbits. At least an order of magnitude more gas particles are needed to significantly improve this situation. While there’s no simple way to anticipate the results of such ambitious calculations, it’s plausible that disk-like structures like those reported in this paper will continue to appear in future experiments.

## 3 ENCOUNTERS

Direct and retrograde versions of a close passage are compared in Figure 1. Encounter DIR 1:1 C (top) resembles direct encounters described in other studies (Barnes & Hernquist 1991, 1996; Mihos & Hernquist 1994, 1996). Although the disks interpenetrate, only a modest fraction of the gas actually collides with gas from the other disk. The stellar and gaseous components both respond to the tidal forces by developing extended bridges and tails (Toomre & Toomre 1972). In the aftermath of such passages, much of the gas is rapidly driven into the central regions; these inflows are driven by gravitational torque between the stellar and gaseous bars formed in tidally perturbed disks (Combes, Dupraz, & Gerin 1990; Barnes & Hernquist 1991; Mihos & Hernquist 1996). Encounter RET 1:1 C (bottom) illustrates a violent hydrodynamic interaction. The geometry of this collision insures that a large fraction of the gas suffers strong shocks as the galaxies intersect. By time $`t=1.25`$ (middle), much of the gas in these disks has been swept toward their centers. These inflows are not driven by gravitational torques – the stellar disks are far less perturbed than the gas, and incapable of exerting strong gravitational torques. Instead, they are driven by hydrodynamic forces; the gas loses its spin angular momentum by colliding directly with gas in the other galaxy. Gas which escapes being swept in forms a plume connecting the two galaxies; such structures may have been observed in some deeply interpenetrating encounters (e.g. Condon et al. 1993; Tsuchiya, Korchagin, & Wada 1998).

To illustrate the range of evolutionary histories in this sample of encounters, Figure 2 concisely summarizes two rather different experiments. In these plots, $`E_{\mathrm{rad}}`$ is the total energy lost to radiative processes, while $`\rho `$ is the gas density. As one might guess from the discussion above, encounter RET 1:1 C is the most dissipative of those studied here; net radiative losses amount to $`20`$ percent of the initial binding energy of the entire system. These losses occur in large-scale shocks as the two galaxies plow into each other at $`t1`$, fall back together at $`t3.7`$, and merge at $`t4.3`$. Gas densities increase with each burst of dissipation; by the end of the simulation $`90`$ percent of the gas lies in a barely-resolved disk at the center of the merger remnant. In contrast, encounter POL 3:1 is the least dissipative, losing only $`5`$ percent of its initial binding energy. The first passage at $`t1`$, while close enough to produce definite tidal features, barely registers in the traces of $`E_{\mathrm{rad}}`$ and $`\rho `$. The next passage, at $`t5.0`$, is more dramatic, and the final merger at $`t6.3`$ drives about half of the gas into the central regions; the rest of the gas eventually settles into a warped disk.

Disk formation need not wait until the merger process is completed; reaccretion from tidal bridges and tails can feed high angular momentum gas back into galaxies after any reasonably direct passage. An example is illustrated in Figure 3. Here the upper row shows gas in the larger disk of encounter DIR 3:1 responding after the direct passage of its lighter companion; note the pronounced bar typically formed in such passages. The lower row shows where shocks occur by rendering each particle with intensity proportional to the local dissipation rate $`\dot{u}`$. While the central bar is prominent in the lower images, shocks also develop in the disk surrounding the bar. Some of these shocks form where gas falling back from the tidal tail and bridge impinges on the disk; others occur within the disk, perhaps in response to forcing by the central bar. Between times $`t=1.875`$ and $`t=3.0`$ the reaccreted gas nearly doubles the size of the disk.

Reaccretion of gas from tidal tails is evident in a recent HI study of NGC 4038/9 (Hibbard et al. 2001). Relative to the systemic velocity, most of the tail associated with the northern disk (NGC 4038) is moving away from us, but the gas at the base of the tail has the opposite sense of motion. This indicates that the material at the base has already attained apocenter and is now falling back onto the disk. This reaccreted gas may be fueling the ring of star formation in the disk of NGC 4038.

One consequence of this disk rebuilding is that the gas and stars may be kinematically segregated in tidally disturbed disks. Reaccreted gas seeks out closed orbits, usually close to the disk plane, while reaccreted stars remain on highly eccentric orbits. This can produce dramatic effects in further tidal interactions, possibly explaining some cases in which stellar and HI tails fail to coincide (Hibbard, Vacca, & Yun 2000, Mihos 2001). Figure 4 illustrates the behavior of gas and stars in the second passage of encounter POL 1:1 C. The first frame shows two fairly relaxed gas disks (grey-scale), each embedded in a distended stellar disk (contours). In the second frame the disks have moved past each other, and a pronounced tail extends from the more direct member of the pair. By the third frame the gas and stellar distributions are quite different – the one is no longer a good predictor of the other. For example, the gas in the lower tail crosses several contours of the stellar distribution; evidently the kinematically cold gas has produced a narrow tidal tail, while the warmer stellar component has raised a broader feature. This last frame also shows that the two galaxies are at the point of merging.

## 4 REMNANTS

In broad outline, the mergers resulting from these twelve encounters are similar to those described in earlier studies (Negroponte & White 1983; Barnes & Hernquist 1991, 1996). Table 2 lists merger times $`t_{\mathrm{merger}}`$ for all twelve encounters. Gas usually has little effect on the overall dynamics of orbit decay and merging; of the four 1:1 C encounters, three merge within times $`\mathrm{\Delta }t_{\mathrm{merger}}\pm 0.2`$ of the corresponding stellar-dynamical versions (Barnes 1998). (The exception is RET 1:1 C, where the mechanical friction of the gas hastens the merger by $`\mathrm{\Delta }t_{\mathrm{merger}}0.5`$.) As mentioned in the introduction, gas which suffers strong shocks tends to fall into galactic nuclei before and during the merger, where it creates central concentrations of highly dissipated gas. On the other hand, the gas which is not involved in strong shocks retains a large fraction of its initial angular momentum, and is available to build extended disks like those shown in Figure 5.

Ideally, simulations would be run well past $`t_{\mathrm{merger}}`$. But as already noted, the later stages of these calculations are quite slow. Practical considerations forced me to stop at the times $`t_{\mathrm{end}}`$ listed in Table 2. Most remnants had $`t_{\mathrm{end}}t_{\mathrm{merger}}2\pm 1`$ time units to relax; scaling the initial disks to the Milky Way, the galaxies in these experiments merged $`10^8`$ to $`10^9\mathrm{yr}`$ before the end of the calculations. These remnants may thus be compared to late-stage Toomre-sequence (Toomre & Toomre 1972; Toomre 1977) galaxies or to IR-luminous merger remnants (Sanders & Mirabel 1996).

### 4.1 Structure

Figure 6 summarizes the distribution and kinematics of gas at the end of each simulation. These plots show logarithmic radii $`\mathrm{log}r`$ vs. radial velocities $`v_r`$ of gas particles. The distributions are roughly symmetric with respect to the horizontal axis $`v_r=0`$ at small radii, but distinctly asymmetric at larger radii. The narrow diagonal streaks seen at large radii represent tails extracted from the disks before they merged. Almost every pericentric passage yields a pair of tails, but the amount of tail material and its subsequent evolution depend on the geometry of the passage. Thus at the instants shown in Figure 6 the remnant produced by RET 3:1 has almost no gaseous tail material, while remnant DIR 3:1 still shows three distinct tails (and the remains of a fourth), each originating from a different disk at a different passage. The dashed curves in the upper right of each plot represent the nominal escape velocity $`v_{\mathrm{esc}}\sqrt{2G(M_1+M_2)/r}`$. Some tail material has enough velocity to escape, but almost all falls back eventually. Tail material with $`v_r<0`$ has passed apogalacticon and is now falling back; as this material reenters the remnant, the gas encounters shocks and begins evolving toward closed orbits.

If the infalling gas has relatively little angular momentum it may plunge deep into the remnant and directly assimilate into a disk. In many cases, however, the tail gas has so much angular momentum that it initially side-swipes the central disk and subsequently follows eccentric, looping paths before falling back yet again; a vivid example is provided on the right in Figure 5. The gas on these looping orbits populates transition zones between asymmetric and symmetric $`v_r`$ distributions in Figure 6. These zones are labeled by the letter “L”; the various POLar encounters all produce clear examples of this behavior. A given gas particle may settle into a closed orbit after executing a few loops, but the loop structures themselves are replenished by continued infall.

At radii where the gas distribution is nearly symmetric with respect to $`v_r=0`$ lie the gas disks formed by this infall process. Gas on nearly-circular orbits has $`v_r0`$, while gas on elongated orbits has a large $`v_r`$ dispersion. Thus in Figure 6 the vertical spread of the points at a given radius gives information on the shapes of gas orbits. For example, remnants such as POL 1:1 C, INC 1:1 C, and INC 1:1 D have large amounts of gas with fairly narrow $`v_r`$ distributions; this gas lies in rotating disks, each labeled “D” in the plots. At smaller radii, remnants DIR 1:1 D, DIR 3:1, POL 3:1, and INC 3:1 have zones, labeled “B” in the plots, with large, symmetric spreads in $`v_r`$. These arise where central triaxial or bar-like structures dominate the potential and the gas consequently moves on elongated orbits.

All twelve remnants have nuclear gas clouds characterized by very high densities, often of order $`\rho 10^6`$ model units. These nuclei, which are labeled “N” in Figure 6, usually have small radial velocity dispersions; most are supported by a combination of pressure and circular motion. But remnants DIR 3:1, RET 3:1, POL 1:1 C, and INC 1:1 C are different; in these cases the nuclear gas lies in rotating elliptical disks. These nuclei are barely resolved in the present simulations, so their detailed dynamics are not well-determined. However, the total masses and angular momenta of nuclear clouds should be more reliable since these are determined by events on somewhat larger scales. Several of these nuclear gas disks are kinematically decoupled with respect to the remnants they inhabit.

The gas fractions in the nuclei, bars, disks, loops, and tails of all remnants are listed in Table 2. Here and in what follows, each of these components is simply defined as all gas within an appropriate range of radii measured from the center of each remnant. These radii often correspond to obvious features in Figure 6; for example, the “D” (disk) component in remnant DIR 1:1 C is delimited by the abrupt end of the infalling tails at $`\mathrm{log}r0.33`$ on the one hand, and by a marked change in the density of points at $`\mathrm{log}r1.92`$ on the other hand. In a few cases the divisions between components are less obvious; in particular, it’s hard to define the boundaries between disks and outer loops from the data in Figure 6. To fix these boundaries more accurately, I used an interactive viewing program to inspect the gas distribution in three dimensions; on such a display, as in the right-hand side of Figure 5, the distinction between nearly circular disks and elongated loops is unmistakable.

As Table 2 shows, the initial orientation of the colliding disks has a definite effect on the distribution of gas in the remnants. Encounters in which both disks are tilted with respect to the orbital plane yield remnants with massive disks, while encounters in which a disk lies in the orbital plane drive more gas to the nuclear regions. Nuclei contain an average of $`85`$ percent of the gas in the remnants of DIRect and RETrograde encounters; the most massive nuclei, found in the RETrograde encounters, are produced by violent hydrodynamic forces as illustrated in Figure 1. In contrast, only $`63`$ percent of the gas lies in the nuclei of the POLar and INClined remnants; the balance is found in the extended bars, disks, and loops surrounding these nuclei.

Figure 7 displays rotation directions for the gas in the remnants. In each plot, an area-preserving map is used to transform the unit sphere $`𝐒`$ of all possible spin directions onto the unit circle $`𝐂`$; the direction of the initial orbit’s angular momentum vector $`𝐉_{\mathrm{orbit}}`$ always maps to the origin of $`𝐂`$. Any area-preserving map from the sphere to the plane must create some distortion; here, the direction exactly opposite to $`𝐉_{\mathrm{orbit}}`$ maps to the entire edge of the circle, and nearby regions are stretched into elongated crescents. However, the area-preserving property of this map implies that an isotropic ensemble of directions is mapped to a uniform distribution within $`𝐂`$. The images represent the angular momenta directions of gas particles; each particle’s angular momentum vector, measured with respect to the center of the remnant as defined by the peak gas density, was projected onto $`𝐒`$ and thence mapped to $`𝐂`$, and the resulting distributions were smoothed with an iterated boxcar. In these images, color is used to indicate radii, with grey for the nuclei, blue for small radii, green and yellow for intermediate radii, and red for large radii.

As Figure 7 shows, the remnants have diverse kinematics. Many of the structures seen here can be identified with components delineated in Figure 6. Gas in the nuclei, shown in grey, produces the broadest distributions, but each nucleus has a well-defined net spin direction which is marked by a “N” in these plots. Likewise, the “B”, “D”, and “L” symbols mark the net spin directions of the bars, disks, and loops; these coincide with the angular momenta of gas particles, shown in blue and green. On the other hand, the spin vectors of the initial disks, labeled “1” and “2”, generally don’t correspond to favored spin directions in the remnants, but often appear associated with material in tidal tails, here shown in yellow and red.

A few remnants display fairly simple structures. For example, remnants of DIRect encounters tend to have nuclei and disks with reasonably well-aligned rotation, basically because the spins and orbital motion of the initial encounter reinforce each another. But in many cases the gas has some kind of kinematic misalignment. All remnants of RETrograde passages have disjoint angular momenta distributions, with counter-rotating nuclei and disks. This is not surprising given the rather contrived initial conditions of these mergers; the distribution of gas between direct and retrograde rotation reflects the competition of spin and orbital angular momenta in these encounters. Several other remnants have gas nuclei with dramatic kinematic misalignments: remnant POL 1:1 C has a counter-rotating nucleus, while remnants INC 1:1 C and D have nuclei which rotate about axes roughly perpendicular to the axes of their disks.

As noted above, many remnant disks have warps. A disk which warps progressively along a single direction appears as a roughly linear feature with a smooth color gradient in Figure 7, as in the plot for remnant INC 1:1 C. In other disks the warp’s direction changes with radius, perhaps as a result of differential precession; such warps produce curved forms as in the plot for remnant POL 1:1 C, or even circular features like the one in RET 1:1 D. Most of these disks are fairly flat at small radii, as if the inner disks are locked into a single plane by self-gravity. Warps in the outer disks of these remnants may be excited by the return of tail material with misaligned angular momentum. In Figure 7 these remnants exhibit fairly smooth transitions from the disks (blue or green) to the tails (yellow or red); thus the angular momentum of the gas changes in a continuous manner with radius. The tails contain large amounts of angular momentum; as this material falls back it could exert strong torques, creating long-lived warps in remnant disks.

Remnant DIR 1:1 C provides the example of a precessing inner disk shown in Figure 8. This disk, bounded by the dotted line at $`\mathrm{log}r1.2`$ in the upper-left panel of Figure 6, begins forming as soon as the nuclei of its progenitors have merged and the gravitational potential has settled down. By $`t=4.5`$ it’s seen roughly edge-on extending from eight to two o’clock; it has precessed counter-clockwise by $`90^{}`$ by $`t=6`$. During this time the inner disk remains relatively well-aligned with the nucleus. But between the nucleus and the inner disk lies a ring of material, inclined by $`130^{}`$, which precesses clockwise by $`120^{}`$ between $`t=4.5`$ and $`6`$. If the nucleus and inner disk are actually locked in alignment, they must be coupled gravitationally, since the ring lies between them.

To rough approximation, both the inner disk and extra-nuclear ring in this remnant behave like inclined rotators precessing in an oblate potential. But the middle row in Figure 8 provides evidence that viscous forces also influence the dynamics of the inner disk. These images show the gas particles weighted by their dissipation rate $`\dot{u}`$; the inner disk is quite prominent, as is the gas just beyond the inner disk. The latter exhibits a spiral dissipation pattern, implying that angular momentum is being transported radially. One consequence of this viscous coupling may be the gradual decrease in the inclination of the inner disk.

The inner disk of DIR 1:1 C acquired its tilt in a fairly straightforward manner. Some $`74`$ percent of the gas particles in this component came from the $`i=71^{}`$ galaxy, and the angular momentum of these particles determines the initial orientation of the disk. Some misaligned components in other remnants may be explained in a similar way; in these cases the misaligned material comes preferentially from one progenitor or the other. But the relationship between the originating disk and the final spin direction is not always simple – gravitational and hydrodynamical torques play a large role in determining the final rotation direction of the gas.

### 4.2 Infall

The simulations cover only the early evolution of these merger remnants; they include the violent dissipative events which form the inner parts of gas disks but end shortly thereafter. On longer time-scales, much of the gas lingering in tidal tails will fall back (Hibbard & Mihos 1995), evolve through loop-like structures, and finally seek out closed orbits. It’s computationally expensive to follow this process much past the stage shown here, but counting the bound tail gas, it appears that the DIRect and RETrograde remnants have enough material to form disks containing up to $`25`$ percent of the total gas, while the POLar and INClined mergers can form disks containing $`25`$ to nearly $`60`$ percent.

The time-scale for this infall is easily estimated in the limit where the tail gas is just barely bound to the remnant. In this limit the tail material moves on nearly radial orbits in an approximately Keplerian potential. Suppose that a tail is formed during a pericentric passage at time $`t_{\mathrm{peri}}`$. A parcel of tail material with specific binding energy $`E<0`$ reaches apogalacticon at radius $`r_{\mathrm{apo}}G(M_1+M_2)/|E|`$, and falls back at time

$$t_{\mathrm{ret}}t_{\mathrm{peri}}+\frac{\pi G(M_1+M_2)}{\sqrt{2}|E|^{3/2}}.$$
(2)
In this model there’s a one-to-one relationship between the binding energy distribution of the tail material and the rate at which this material is reaccreted. Let $`dm/dE`$ be the differential distribution of binding energy of tail gas; then the reaccretion rate is

$$\frac{dm}{dt}=\frac{dm}{dE}\frac{dE}{dt_{\mathrm{ret}}}\frac{dm}{dE}(tt_{\mathrm{peri}})^{5/3}.$$
(3)
As Figure 6 shows, some tails are entirely bound; these fall back completely after a finite time. On the other hand, a tail which extends to $`E>0`$ generally has a distribution $`dm/dE`$ which is approximately constant for $`E0`$. This follows because the tidal interaction has no way to “pick out” the escape energy; thus $`dm/dE`$ will be nonzero and continuous near $`E=0`$ (e.g. White 1987; Jaffe 1987). Thus if some tail gas has enough energy to escape then the reaccretion rate $`dm/dtt^{5/3}`$ as $`t\mathrm{}`$.

This model gives a crude but serviceable description of the infall of tail material in remnant DIR 1:1 C. Figure 9 shows infall times for particles from the tails created at the first passage ($`t_{\mathrm{peri}}=1`$); the horizontal axis is the measured time of infall, while the vertical axis is the time predicted by (2) using specific binding energies determined at time $`t=4`$. This comparison was made using collisionless particles from the stellar disks: gas is deflected by shocks, so its exact instant of perigalacticon is harder to define. But both gas and stars travel on nearly ballistic trajectories while in tails, so the rough agreement seen here should also hold for the gas. The scatter between the predicted and measured infall times probably arises because actual orbits are not radial and the true potential is not Keplerian. At later times the simple analytic model should grow more accurate since late-returning material falls back from ever-larger distances. Figure 10 shows the infall rate for gas from the first-passage tails of DIR 1:1 C, predicted using binding energies measured at $`t=4`$. Because of the way this sample of tail particles was defined, the infall rate is depressed for $`\mathrm{log}t0.7`$. At later times, however, the infall rate evolves almost precisely as $`t^{5/3}`$. This is consistent with the distribution of binding energies $`dm/dE`$ of gas in these tails, which varies by less than $`\pm 20`$ percent for $`1.25<E<0.125`$.

If the infall time of tail gas is basically determined by its specific binding energy $`E`$, its evolution once it falls back into the remnant depends largely on its specific angular momentum $`J`$. Studies of disk formation (e.g. Mestel 1963; Fall & Efstathiou 1980) often assume that each parcel of disk material conserves angular momentum during proto-galactic collapse. This can be tested in the simulations; Figure 11 compares specific angular momenta of gas particles at $`t=4t_{\mathrm{merger}}`$ and $`t=6=t_{\mathrm{end}}`$ in remnant DIR 1:1 C. This plot shows two main components: a broad spray of points in the lower left, and a more or less diagonal distribution extending to the upper right. The former is the nucleus, while the latter include the extended disk and tidal tails. The nucleus captures some gas between $`t_{\mathrm{merger}}`$ and $`t_{\mathrm{end}}`$, but most of the gas which has enough angular momentum at the time of the merger – say $`\mathrm{log}J1.4`$ – remains outside the nucleus and approximately conserves the magnitude of its specific angular momentum.

Detailed conservation of angular momentum can be used to predict the radial distribution of the gas at later times. For simplicity, assume that the gas eventually finds circular orbits in the potential $`\mathrm{\Phi }(r)`$ of a spherical mass distribution. Then a gas particle’s specific angular momentum $`J`$ and orbital radius $`r`$ are related by

$$J=\sqrt{r^3d\mathrm{\Phi }/dr}.$$
(4)
Let $`dm/dJ`$ be the differential angular momentum distribution of the gas; then the predicted radial distribution in the final disk is

$$\frac{dm}{dr}=\frac{dm}{dJ}\frac{dJ}{dr}.$$
(5)
To predict how the disk’s radial profile evolves as infall proceeds, (2) may be used to restrict $`dm/dJ`$ to only the gas which has fallen in before a given time.

Figure 12 illustrates the prediction of radial coordinates for gas particles in remnant DIR 1:1 C using (4) and the model potential in Appendix A. Three different scatter-plots are presented; in every one the vertical axis shows radii of gas particles which lie in the range $`1.6<\mathrm{log}r<0.33`$ at time $`t=6`$. In the top plot the horizontal axis shows radii for these particles calculated from (4) using angular momenta $`J`$ measured at $`t=6`$; the points fall almost precisely on the diagonal line, confirming that the gas particles are on nearly circular orbits and that the mass model is adequate. In the middle plot the same equations are used, but the angular momenta $`J`$ are measured at $`t=4`$, which is just after the merger. If individual gas particles conserved $`J`$ precisely, this plot would be identical to the top one; the actual distribution is considerably broader, but still more or less straddles the diagonal line. Finally, in the bottom plot the horizontal axis shows the radii measured just after the merger at $`t=4`$; this distribution is broader still and distinctly offset from the diagonal line, showing that significant infall has taken place between $`t=4`$ and $`t=6`$.

In sum, Figure 12 represents a qualified success for this simple model of disk formation by infall from tails. Given a model potential $`\mathrm{\Phi }(r)`$, the radial coordinate of a particle of disk gas can be calculated with some accuracy if its final angular momentum is known. If the angular momentum at an earlier time is given instead, it’s still possible to roughly predict the final radius. Such predictions show some scatter, but they’re better than simply taking radii at the moment of merger as estimates of final radii. The scatter could be reduced if the torques acting on the gas between $`t_{\mathrm{merger}}`$ and $`t_{\mathrm{end}}`$ were taken into account – but this is beyond the scope of the present model.

Figure 13 illustrates the use of (4) to predict the future growth of the disk in remnant DIR 1:1 C. Here the dotted curve shows the actual radial distribution of all gas particles with $`E<1`$ at $`t=6`$; this subset includes the disk and some of the infalling tail material. The three solid curves show predicted radial distributions of gas particles with $`E<1`$, $`E<0.5`$, and $`E<0`$. At radii $`\mathrm{log}r0.4`$ all particles are on essentially circular orbits and all four curves agree. At larger radii the actual distribution is more extended than the predicted distribution for $`E<1`$ since outlying gas has not yet had a chance to settle into circular orbits. The distribution predicted for $`E<1`$ has a fairly sharp break at $`\mathrm{log}r0.2`$ and essentially vanishes at $`\mathrm{log}r0`$. This reflects the fact that there is a fairly well-defined ceiling on the angular momenta of gas particles with $`E<1`$; binding energies and angular momenta in tails are correlated, with the less tightly bound tail material having more angular momentum. Gas with $`1<E<0.5`$ has more angular momentum than gas with $`E<1`$, and therefore settles at larger radii; gas with $`0.5<E<0`$ settles at larger radii still. Thus the disk will grow from the inside-out, with reaccretion constantly depositing gas at the outer edge of the disk. This can be seen in Figure 8; the top row of images show the radius of the disk increasing with time, and the middle row shows that the periphery of the disk has a high rate of dissipation as infalling material collides with disk gas.

## 5 DISCUSSION

Gas disks like those described here were briefly reported in earlier simulations of dissipative mergers between disk galaxies (Hernquist & Barnes 1991; Barnes & Hernquist 1996). These studies used the same general methodology as the present work, but differed in several specifics. First of all, in earlier work the thermal evolution of the gas was explicitly computed, including radiative cooling with a cut-off at $`T=10^4\mathrm{K}`$. This had little consequence for the gas in regions of moderate to high density where cooling times are short; however, the gas tails in the present experiments are broader and smoother since they are not cooled by adiabatic expansion. Second, the earlier studies used King (1966) distributions for the bulges and halos of the initial galaxy models. As a result, these models may have been more susceptible to bar instabilities (Mihos & Hernquist 1996), and their relative orbits decayed appreciably faster (Barnes 1998). Third, the earlier simulations used the TREESPH code (Hernquist & Katz 1989); the code used in the new simulations borrows several features from TREESPH, but it adopts a different smoothing prescription for gravitational forces, a more efficient and robust tree algorithm, a fixed number of neighbors in SPH summations, and a different time-step algorithm. Despite using somewhat different models and codes, these studies all yield broadly similar results; the process of disk formation appears to be insensitive to such details.

Simulations of polar ring formation in galaxy mergers have produced remnants with very large rings or annuli (Bekki 1998). These rings arise from a relatively restricted set of initial conditions which don’t overlap with the ones considered here.

Disk formation is also seen in simulations of hierarchical clustering which include gas as well as cold dark matter (e.g. Evrard, Summers, & Davis 1994; Navarro & Steinmetz 1997; Domínguez-Tenreiro, Tissera, & Sáiz 1998). These disks are generally identified with spiral galaxies, and they form largely from pristine gas which was not previously part of another galaxy. In common with the disks seen in the present experiments, disks formed by hierarchical clustering generally grow from the inside out.

### 5.1 Origins of counter-rotation

A rather striking result from this study is the large number of remnants with counter-rotating or otherwise misaligned nuclei. In some cases, counter-rotation is easily explained; for example, it’s no great surprise that the RETrograde encounters all produce remnants with counter-rotating components. Likewise, some misaligned components appear to contain gas originating almost exclusively from one of the two progenitors, and it’s likely that such “unblended” structures retain some memory of the orientation of their parent disks. But the origin of the misalignment is not always so evident.

The nuclear disk in remnant POL 1:1 C is a case in point. In terms of overall structure this remnant is very similar to one previously described (Hernquist & Barnes 1991); in both a fairly massive nuclear gas disk counter-rotates with respect to an outer gas disk and stellar component. This resemblance may not seem surprising, since the same initial disk geometry ($`i_1=71^{}`$, $`\omega _1=90^{}`$, $`i_2=109^{}`$, $`\omega _2=90^{}`$) was used in both experiments. But the initial orbit used in the previous study was wider, with a first pericenter at $`r_{\mathrm{peri}}=0.4`$. Of the encounters in the present study, the one with initial conditions most closely matching those used previously is POL 1:1 D, and POL 1:1 D does not yield a counter-rotating nucleus. It’s not clear why this encounter doesn’t produce a counter-rotating nucleus like the one seen in the earlier study, while it’s close relative POL 1:1 C provides a very good match. One possibility is that the halos used in the present study are a bit less effective at promoting orbital decay; as a result, the merging nuclei arrive at their final encounter with slightly more orbital angular momentum than they had in the earlier experiment.

The mechanism which forms counter-rotating nucleus in objects like remnant POL 1:1 C remains obscure. The nucleus of this remnant is already counter-rotating when it forms, and it’s difficult to determine when and how this material first acquired its rotation; angular momentum can only be measured with respect to a center, and before the nucleus forms the system has no unique center. This ambiguity in defining the center could actually be part of the puzzle – material co-rotating with respect to the center as defined in one way may be counter-rotating with respect to the center defined in another. Large-scale hydrodynamic interactions seem to be important; many of the counter-rotating or misaligned nuclei in the present sample form when galaxy pairs arrive at their second passage still dressed in extensive gas disks. A better understanding of the roles of gravitational and hydrodynamic torques in the formation of counter-rotating nuclei may require calculations with higher spatial resolution than available at present.

### 5.2 Observational connections

The results presented here invite comparisons with observations of gas disks in merger remnants. As mentioned in the introduction, NGC 7252 is an obvious example; the central disk of ionized and molecular gas in this system is about the same size as the disks produced in these experiments, and the peculiar kinematics of this disk may be due to a strong warp of the kind seen in many of these simulations. In addition, HI velocities confirm that gas is falling back into the remnant from the tidal tails (Hibbard et al. 1994). NGC 7252 has been modeled as the result of a close and fairly direct encounter between two disk galaxies of comparable mass (Hibbard & Mihos 1995). This model nicely reproduced the morphology and kinematics of the tails, but the calculations did not include a dissipative component which could form a central disk. It’s worth repeating this calculation with a combined N-body/SPH code; a single model reproducing both the tails and the central disk would be a fairly impressive accomplishment, and might help constrain the gas content of NGC 7252’s progenitors.

The peculiar elliptical galaxy NGC 3656 (Balcells 1997, Balcells et al. 2001) may be a second example. This galaxy contains shells and a pair of faint tails which suggest a merger of two disk galaxies. It also has a star-forming dust lane which corresponds to an extended disk of HI. This disk is visibly warped, and the outer edge of the disk is kinematically contiguous with gas at larger radii which may be falling in from the tidal tails. The estimated gas accretion rate in NGC 3656 is an order of magnitude lower than in NGC 7252; this may imply that NGC 3656 is at a later stage in its evolution than NGC 7252.

NGC 5128, an elliptical galaxy with an active nucleus and an extended disk of dust and gas, may represent an even later stage in the evolution of gas-rich merger remnants. Although the warped disk in this system has been interpreted as the remains of an accreted gas-rich satellite (e.g. Malin, Quinn, & Graham 1983), there are some grounds to suspect that this galaxy is the result of a fairly major merger (e.g. Schweizer 1998). These include the ripples visible in deep optical images (Malin et al. 1983), the HI fragments surrounding the galaxy (Schiminovich et al. 1994), and the misaligned rotation revealed by kinematics of planetary nebulae (Hui et al. 1995).

Nuclear gas disks and rings seem to be very common in merging galaxies selected by infrared luminosity (e.g. Downes & Solomon 1998). In remnants with single nuclei the gas tends to rotate more rapidly than the stars, indicating some degree of kinematic decoupling (Genzel et al. 2001), though it’s not clear if the extreme decoupling in remnant POL 1:1 C has yet been seen in the observations. Disks in systems which appear to still have two nuclei are not so easily explained. Molecular-line observations of Arp 220 and NGC 6240 have been interpreted in terms of a rotating disk located between the nuclei (Scoville, Yun, & Bryant 1997; Tacconi et al. 1999; Tecza et al. 2000). In the simulations, gas driven inward before the galaxies merge always accumulates in disks or bars around the individual nuclei; moreover, an incipient disk between the nuclei would be torn apart by gravitational fields as the orbits of the nuclei decay. An inter-nuclear disk might be stable if it was more massive than either nucleus (Tecza et al. 2000), but the amount of gas required seems quite extravagant, and earlier reports of a peak in the stellar velocity dispersion between the nuclei of NGC 6240 (Lester & Gaffney 1994; Doyon et al. 1994) are not supported by recent HST observations (Gerssen et al. 2001). On the other hand, the interpretation of Arp 220’s nuclei as bright spots in “a warped molecular gas disk” (Eckart & Downes 2001; Downes & Solomon 1998) seems entirely consistent with the present numerical results.

### 5.3 Galaxy transformation

Simulations indicate that extended gas disks can form in mergers of spiral galaxies. These disks contain between $`20`$ and $`60`$ percent of the total gas in the initial galaxies, and may extend to several times the remnant half-light radii. It’s worth noting that the initial galaxy models used in these experiments were not particularly gas-rich, and that the gas started with the same distribution as the disk stars. In real spiral galaxies the atomic gas is more extended than the stars; mergers of such systems should form remnants with even larger and more massive gas disks. Moreover, the gas content of disk galaxies is generally expected to increase with redshift; remnants with disks containing $`10`$ to $`15`$ percent of their luminous mass would very likely result by doubling the gas fraction in the initial galaxy models.

These disks, if subsequently converted to stars, would be fairly hard to detect photometrically (Rix & White 1990) unless viewed from a favorable orientation (e.g. Scorza & Bender 1990). On the other hand, disks could significantly influence the observed kinematics of merger remnants. Disks are kinematically cold, so stars formed in such disks would add narrow features to the stellar absorption-line spectra of early-type galaxies. Moreover, disks rotate faster than pressure-supported components, so these narrow features will be systematically offset from the broader profiles due to the rest of the galaxy. The combined velocity profile of a disk plus spheroid thus appears asymmetric, with a steep prograde wing and a shallow retrograde wing (e.g. Franx & Illingworth 1988; Bender 1990; Rix & White 1992; van der Marel & Franx 1993). Such profiles actually appear to be fairly typical in elliptical galaxies with measurable rotation (Bender, Saglia, & Gerhard 1994).

The observed velocity profiles may be compared with those predicted for remnants of purely dissipationless mergers between disk galaxies (Bendo & Barnes 2000, Naab & Burkert 2001). As a rule, simulated equal-mass mergers produce remnants with the opposite velocity profile asymmetry – that is, they have shallow prograde wings and steep retrograde wings. For unequal-mass mergers the situation is not so clear – one study finds profiles similar to those produced in equal-mass mergers, while the other reports cases which exhibit the same sense of asymmetry seen in observations of early-type galaxies. Naab & Burkert (2001) suggested that the formation of extended disks is necessary if mergers are to account for the observed kinematics of elliptical galaxies, but worried that the gas might not retain enough angular momentum to build extended disks. Given the contradictory results on unequal-mass mergers, it may be premature to insist that all mergers must form disks. On the other hand, the present work indicates that disks of the requisite size can form quite easily in both equal-mass and unequal-mass mergers.

If the above explanation for line profiles in elliptical galaxies is correct, one might expect different spectral lines to show different kinematic signatures. The lines produced by the younger stars making up a disk should be narrower than, and systematically offset from, the lines produced by the older stars in the spheroid. This might be tested with high signal-to-noise spectra obtained using large optical telescopes.

## ACKNOWLEDGMENTS

I thank the referee, A. Burkert, for a constructive report, and M. Balcells, L. Hernquist, J. Hibbard, D. Sanders, and F. Schweizer for valuable discussions. This work was supported in part by STScI grant GO-06430.03-95A.

## Appendix A REMNANT POTENTIAL

A simple model of the potential may be constructed using a Plummer model to represent the nucleus and a Hernquist model to represent the rest of the remnant:

$$\mathrm{\Phi }(r)=\frac{GM_{\mathrm{nuc}}}{(r^2+a_{\mathrm{nuc}}^2)^{1/2}}\frac{GM_{\mathrm{rem}}}{r+a_{\mathrm{rem}}}.$$
(6)
For remnant DIR 1:1 C, I took $`M_{\mathrm{nuc}}=0.034`$ to be the mass within $`r<0.005`$, $`a_{\mathrm{nuc}}=ϵ=0.0125`$ to be the “softened” radius of the nucleus, $`M_{\mathrm{rem}}=M_1+M_2M_{\mathrm{nuc}}=2.466`$ to be the mass of the rest of the remnant, and $`a_{\mathrm{rem}}=0.2128`$ to be the radius enclosing a quarter of the remnant’s mass exclusive of the nucleus.


text:
# Sensitivity to Initial Conditions in Quantum Dynamics: an Analytical Semiclassical Expansion

## The zeroth order approximation: The mean field approximation (MFA) and the classical limit.

Our zeroth order approximation is defined in the following way:

a) We consider here states of the form of products of coherent states and a phase

$$|\psi (t)=\mathrm{exp}(i\eta (t))|x(t)|y(t),$$
(2)
where the minimum uncertainty coherentes states $`|z(t)`$ are defined by $`|z(t)=𝒟(z(t))|0`$ with

$$𝒟(z(t))=\{\begin{array}{c}\mathrm{exp}[z(t)a^{}z^{}(t)a]\hfill \\ \mathrm{exp}[\frac{\mathrm{atan}|z(t)|}{|z(t)|}(z(t)J_+z^{}(t)J_{})],\hfill \end{array}$$
(3)
for h(3) and su(2) respectively, and the fiducial state $`|0`$ is the Fock state $`|n=0`$ for h(3) and the $`J_z`$ eigenstate $`|J,J`$ for su(2).

b) Linear combinations are not allowed as initial conditions (this circumvents problems with the superposition principle). Also this MFA for the systems (1) will leave invariant the manifold of the chosen set of trial functions.

At this point it is important to mention that the algebra h(3) could be easily enlarge to include of harmonic oscillator algebra with little effort, but the extra terms give rise to nonlocality effects (squeezing dynamics) already at this lowest order which we would like to avoid. Notice also that generalization to n degrees of freedom m-linear hamiltonians ($`mn`$) is straightforward.

Now we solve the mean field Schrödinger equation

$$(i_t\overline{H})|\psi (t),$$
(4)
where $`\mathrm{}=1`$ and the mean field hamiltonian (MFH) is

$$\overline{H}=1_AH_B(x(t))+H_A(y(t))1_B,$$

with $`H_B(x(t))=x(t)|H|x(t)`$. Using the explicit form of the hamiltonian (1) we find

$$H_A(t)=\underset{i}{}a_i(t)A_i+f_A(t),$$
(5)
where $`a_i(t)=\alpha _i+_j\gamma _{ij}y|B_j|y`$ and $`f_A(t)=_j\beta _jy|B_j|y`$. The expression for $`H_B(t)`$ is completely analogous to the expression for $`H_A(t)`$. Since the terms $`f_A(t)`$, $`f_B(t)`$ give rise to global phases we neglect them. Due to the structure of the MFH the solutions of eq. (4) will preserve the form (2), and their time dependence will be completely specified by the solutions of the following equations of motion

$`\dot{z}`$ $`=`$ $`ia_0zia_+(\text{h(3)})`$ (6)
$`\dot{z}`$ $`=`$ $`ia_+ia_0z+i(a_+)^{}z^2(\text{su(2)})`$ (7)
$`\dot{\eta }_z`$ $`=`$ $`0|𝒟^{}(z(t))\left(i_th(t)\right)𝒟(z(t))|0.`$ (8)
Observe that the total phase $`\eta (t)`$ is the sum of the partial phases $`\eta _x(t)`$ and $`\eta _y(t)`$. The nonlinearity of these equations arise from the self consistency of MFA. If the labels are scaled as $`z=\sqrt{4J}Z`$ for h(3), $`z=Z`$ for su(2), and time as $`t=t_c/(4J)`$, then the equations for $`Z(t_c)`$ will become independent of $`J`$ and correspond to the classical limit of Heisenberg’s equations of motion. From equation (8) it follows that the phases $`\eta _z(t)`$ are generalized actions of the coherent state $`|z(t)`$. They are of course absent from the classical limit, but will be crucial for the quantum corrections.

Thus, we have shown that systems with hamiltonian (1), in the mean field approximation satisfy all the requirements we wanted: labels with classical physical meaning, no superposition principle, minimum quantum nonlocality effects, and hamiltonian equations for labels which coincide with both the exact Schrödinger equation and the Heisenberg equations of motion.

## Quantum corrections.

Now we turn to the corrections to the MFA. The exact Schrödinger equation for the whole system (with the tilde indicating the Schrödinger picture)

$$i_t\stackrel{~}{|\psi (t)}=(\overline{H}(t)+\stackrel{~}{\mathrm{\Delta }}(t))\stackrel{~}{|\psi (t)},$$

where $`\stackrel{~}{\mathrm{\Delta }}(t)=H\overline{H}(t)`$, can be written in the MFA interaction picture (MFAIP) as (the absence of the tilde indicating the MFAIP)

$$i_t|\psi (t)=\mathrm{\Delta }(t)|\psi (t)=\underset{ij}{}\gamma _{ij}(A_i(t)A_i(t))(B_j(t)B_j(t))|\psi (t),$$
(9)
again up to time dependent terms that give rise to a global phase. Here we have defined $`A(t)`$ and $`A(t)`$ as

$$\begin{array}{cc}\hfill A(t)& =\overline{U}^{}(t,0)A\overline{U}(t,0),\hfill \\ \hfill A(t)& =0|𝒟^{}(x(t))A𝒟(x(t))|0,\hfill \end{array}$$

where $`\overline{U}(t,0)`$ is the evolution operator of the MFA, the product of the evolution operators for each degree of freedom. The equations for $`B(t)`$ and $`B(t)`$ are completely analogous. Equation (9) possesses a natural expansion

$$|\psi _t=(1i_0^t𝑑t_1\mathrm{\Delta }_{t_1}_0^t𝑑t_1_0^{t_1}𝑑t_2\mathrm{\Delta }_{t_1}\mathrm{\Delta }_{t_2}+\mathrm{})|I,$$
(10)
where $`|I`$ is the initial state $`|x(0),y(0)`$. Now, let us see that all the terms in this expansion are readily written in terms of classical quantities. For example, the first correction term gives

$$|\psi (t)=|IiC(t)|D$$

where C(t) is a function to be defined below and $`|D`$ is a generalized coherent state, orthogonal to the initial coherent state, given by $`|D=|D(z(0))=𝒟(z(0))|1`$ with the reference state $`|1`$ given by the Fock state $`|n=1`$ for h(3) and by $`|J,J+1`$ for su(2). We immediately see that the first correction already introduces all of the effects we avoided at the mean field level: superposition of states, nonlocality and entanglement. State $`|D`$ is sometimes called a doorway state. In order to get the first order result we made use of the following identity which is a result of the present group(s) structure. In fact, it is this relation which enables one to express all order corrections in terms of classical trajectories and generalized coherent states.

$$A_iD(x)=D(x)\left(\underset{k}{}g_{ik}^A(x)A_k+k_i^A(x)\right),$$
(11)
where $`g_{ik}(x)`$ are functions of $`x`$ specific to each one of the groups in question. Similar relations hold for the degree of freedom B. The coefficient $`C(t)=_0^tc(t_1)𝑑t_1`$, where

$$c(t)=\sigma e^{i(S_0S_1)(t_1,0)}\underset{ij}{}\gamma _{ij}g_{i+}^A(x(t_1))g_{j+}^B(y(t_1))$$

is clearly only a function of classical trajectories and the corresponding actions, $`S_0(t,0)=_0^t𝑑t_1(\eta _x(t)+\eta _y(t))`$ and $`S_1(t)=S_1[x(t)]+S_1[y(t)]`$ with

$$S_1[z(t)]=_0^t𝑑\tau 1|𝒟^{}(z(\tau ))\left(i_\tau h(\tau )\right)𝒟(z(\tau ))|1.$$

It is clear that the remaining corrections can also be written in terms of classical trajectories (and actions), but their quantum content will not be as transparent as in the leading correction. In fact, the second correction can be written as having a term proportional to the initial coherent state, but also terms proportional to generalized coherent states $`|x_0,y_1`$, $`|x_1,y_0`$, and $`|x_2,y_2`$, where the subindices refer to the fiducial states.

## Sensitivity to Initial Conditions: A Formal Nonperturbative Result

Classically one of the basic ingredients to define chaos is the high sensitivity to initial conditions. A formalization of this condition is heavily based on the notion of distance between trajectories. In establishing a quantum counterpart of this condition it is important to introduce a quantum measure of distance between states. A natural measure is given by the square modulus of the scalar product.222Some proposals have been made related to scalar product of wavefunction evolved from different hamiltonians (not different wavefunctions). See A. Peres, Phys. Rev. A 30, 1610 (1984); R.A. Jalabert, and H.M. Patawski, Phys. Rev. Lett. 86, 2410 (2001). In what concerns our mean field approximation, the squared modulus of the scalar product between different states of the manifold allows for a direct association of the distance between states with distances between phase space trajectories, since $`|z_1|z_2|^2`$ is given by

$$\{\begin{array}{cc}\mathrm{exp}(|z_1z_2|^2)\hfill & \text{ for h(3),}\hfill \\ (1\frac{|z_1z_2|^2}{(1+|z_1|^2)(1+|z_2|^2)})^{2J}\hfill & \text{ for su(2).}\hfill \end{array}$$
(12)
The important quantum tool which allows us to investigate the sensitivity to initial conditions of the quantum dynamics and eventually make connection to the well known classical limit is the overlap between two time dependent states, which evolved from different initial conditions. It is well known that for unitary evolutions the scalar product is conserved in time. Observe that different initial states correspond to different MFAs, since the MFA is state dependent, due to self consistency. Thus the scalar product between two different initial states is to be written as

$`x^{}(0),y^{}(0)|x(0),y(0)=x^{}(t),y^{}(t)|x(t),y(t)+x^{}(t),y^{}(t)|(\delta _{QC}^{}(t))^{}|x(t),y(t)`$ (13)
$`+x^{}(t),y^{}(t)|\delta _{QC}(t)|x(t),y(t)+x^{}(t),y^{}(t)|(\delta _{QC}^{}(t))^{}\delta _{QC}(t)|x(t),y(t)`$
where we have written $`U(t)`$, the exact quantum time evolution operator as $`\overline{U}^{()}(t)(1+\delta _{QC}^{()})`$, with $`\overline{U}^{()}(t)`$ the MFA evolution operator corresponding to state $`|x(t)^{()},y(t)^{()}`$ and $`1+\delta _{QC}^{()}`$ the evolution operator for the quantum corrections.

The first term on the rhs, $`x^{}(t),y^{}(t)|x(t),y(t)`$ contains the mean field approximation and is given by

$$e^{d(x(t)x^{}(t))/2}e^{i\mathrm{\Phi }(x(t),x^{}(t))}e^{d(y(t)y^{}(t))/2}e^{i\mathrm{\Phi }(y(t),y^{}(t))},$$

where $`\mathrm{\Phi }`$ is some phase which also depends on the classical trajectory and the corresponding group, and functions $`d`$ can be determined by comparison with (12). This matrix element is proportional to the distance between the labels which, in the present case, corresponds precisely to the classical trajectories. Since the exact evolution preserves overlap, the sum of this term with the other three of eq. (13), which contain the quantum corrections, should be conserved in time. Observe that the rate of change of this overlap can have two distinct origins, dephasing and/or change of the modulus. Both changes should be compensated by quantum corrections, but only the later one can be unambigously connected to classical chaos, since, if the system is classically chaotic this distance will exponentially grow and, as a consequence, the overlap involving only the MFA will decrease accordingly. This is a quantum counterpart of the fact that classically chaotic systems exhibit high sensitivity to initial conditions. The corresponding quantum system will exhibit a high sensitivity to the initial state in what concerns the production rate of non unitary quantum corrections to this overlap. We remark that this result is exact and independent of the approximation used. For the argument, however it has been crucial to separate the mean field contribution explicitly. It should also be emphasized that the time scale for the overlap quantum corrections is essentially linked to the Lyapunov exponents. However, other observables will have different time scales, sometimes much shorter, as for example that for the entanglement process. Consequently, entanglement is not always a good measure of classical chaotic behavior, at least for short times, unless the exponential separation of neighbouring classical trajectories also occurs at very short time scales. Should the exponential separation occur at early times, a significant increase in linear entropy will be noticed. In the cases the two time scales are very different this effect, although present, will be rendered less conspicuous by the time development of quantum correlations stemming from the other sources. This will become clear in the example below.

In order to characterize the degree of entanglement we will calculate the idempotency defect (or linear entropy) $`\delta (t)=1\mathrm{Tr}_\mathrm{A}(\rho _A(t))^2`$, where the reduced density $`\rho _A(t)`$ is given by $`\rho _A(t)=\mathrm{Tr}_\mathrm{B}|\psi (t)\psi (t)|`$. Observe that this measure of entanglement does not depend on the picture used to calculate it. Using the expansion (10) up to second order, writing $`|I=|I_A|I_B`$, and calculating the idempotency defect in the MFAIP we obtain to second order

$$\delta (t)=4\mathrm{Re}_0^t𝑑t_1_0^{t_1}𝑑t_2c^{}(t_1)c(t_2).$$

## Application to the classicaly chaotic maser model

The classically chaotic maser hamiltonian

$$H=ϵJ_z+\omega a^{}a+\frac{G}{\sqrt{J}}(a^{}J_{}+aJ_+)+\frac{G^{}}{\sqrt{J}}(a^{}J_++aJ_{})$$

belongs to the class of bilinear hamiltonians (1), where the field (atomic) degree of freedom A (B) is related to the h(3) (su(2)) algebra. Field coherent states are characterized by $`x(t)`$ while spin coherent states by $`y(t)`$. The MFH is of form (5) for each degree of freedom with

$`a_+`$ $`=`$ $`i{\displaystyle \frac{Gy+G^{}y^{}}{\sqrt{J}}},a_{}=(a_+)^{}a_0=\omega ,`$
$`a_+`$ $`=`$ $`i{\displaystyle \frac{Gx+G^{}x^{}}{\sqrt{J}}},a_{}=(a_+)^{}a_0=ϵ,`$
for the field and atomic MFH respectively. Also, the equations of motion for each degree of freedom are equations (6) for the field degree of freedom and (7) for the atomic one. For this model it is a simple matter to give an analytic expression for $`c(t)`$, the key ingredient for evaluation of both the first order correction for the state, and the idempotency defect.

$$c(t)=\frac{\sqrt{2}e^{i(S_{0,J}S_{1,J+1})(t)}}{1+|y(t)|^2}(G^{}Gy^2(t))$$

where $`S_{0,J}`$ ($`S_{1,J=1}`$) is the generalized action for the coherent state with fiducial state $`|n=0|J,J`$ ($`|n=1|J,J+1`$). In this case the doorway state $`|D`$ is given by $`|D=|x_1,y_{J+1}`$.

The time development of the magnitude of the overlap in the MFA between two coherent states centered in the classically chaotic phase space region and also the overlap between other two states chosen in the regular region are shown in fig. (1), and illustrate the effect of the classically chaotic motion on quantum dynamics is dramatic. In effect, for the times the overlap changes appreciably, one can check that the classical trajectories of the chaotic region in question are also correspondingly well set apart. It is clear that for the magnitude of the overlap in MFA the time scale of correlation effects is essentially dictated by classical dynamics. Of course, this needs not hold for other quantum observables. Entanglement, for example, is a quantum property with a smaller time scale. The expected sudden increase in this quantity at the time the modulus of the overlap diminishes, disappears in the midst of contributions of several quantum effects other than the one related to the classical limit (see ref. ). Our analytical approximation for entanglement breaks down for times of the order of 2 for initial conditions both in the chaotic and regular regions. Generalization of these results to other quantum systems is possible, but quantum effects such as nonlocality will be already present at the lowest order, and other analytical approximations should be advanced. Work along these lines is in progress.

The authors are grateful to E.V. Passos and M.P. Pato for fruitful discussions. This work was partly funded by FAPESP, CNPq and PRONEX (Brazil), and Colciencias, DINAIN (Colombia). K.M.F.R. gratefully acknowledges the Instituto de F’ısica, Universidade de São Paulo, for their hospitality and PRONEX for partial support.


text:
# Activating NPPT distillation with an infinitesimal amount of bound entanglement

## I Introduction

The possibility of distillation plays a crucial role in Quantum communication and Quantum information processing (cf. Springer ). Together with quantum error correction it enables all the fascinating applications provided by Quantum information theory in the presence of a noisy and interacting environment. Despite its practical relevance and quite considerable effort in that direction, however, many of the basic questions concerning distillation are yet unanswered. Most notably the question whether a given quantum state is distillable or not, i.e., whether it is possible to obtain pure maximally entangled states from several copies of it by means of local operations and classical communication (LOCC).

A necessary condition for the distillability of a state described by a density matrix $`\rho `$ is the fact that its partial transpose $`\rho ^{T_A}`$, defined with respect to a given product basis by $`ij|\rho ^{T_A}|kl=kj|\rho |il`$, has a negative eigenvalue HHH2 . Except for special cases like states on $`^2^n`$ 2n ; IHa  and Gaussian states Gauss  it is however unclear whether this condition is sufficient as well. There is some evidence presented in IHa ; IBM  that this may not be the case and that there are indeed undistillable states, whose partial transpose is not positive (NPPT). At least there exist $`n`$-undistillable NPPT states for every finite $`n`$, meaning that no LOCC operation on $`n`$ copies leads even to a single entangled two qubit state IHa ; IBM .

However, if we enlarge the class of allowed distillation protocols from LOCC to channels respecting the positivity of the partial transpose, then every NPPT state becomes 1-distillable EVWW ; KLC2  (which can be shown by using entanglement witnesses KLC2 ). Moreover, it is a result from CDK  that these channels can always be stochastically implemented by an LOCC operation where the two parties are given an entangled state with positive partial transpose (PPT) as an additional resource. The latter is known to be bound entangled since the entanglement needed for the preparation of the state cannot be recovered by distillation HHH2 . Nevertheless, PPT bound entangled states can be useful in order to activate the distillability of bipartite NPPT states Terhal ; KLC .

The aim of the present paper is to investigate the limits and requirements of such an activation process. We will show that there exist states with an arbitrary small amount of PPT bound entanglement, which are capable of activating any NPPT state. The required additional resource is therefore universal as well as arbitrary weakly entangled.

## II Preliminaries on symmetric states

One of the key ideas in what follows will be the exploitation of the symmetry properties of states commuting with certain local unitaries. Two well known one-parameter families of such states are the Werner states and Isotropic states, both playing an important role in the sequel.

Werner states W89  acting on a Hilbert space $`=_A_B`$ with dimensions $`dim_A=dim_B=d`$ commute with all unitaries of the form $`UU`$ and can be written as

$$\rho (\alpha )=\left(\mathrm{𝟏}\frac{\alpha }{d}𝔽\right)/(d^2\alpha ),\alpha [d,d],$$
(1)
with $`𝔽`$ being the flip operator defined with respect to some product basis by $`𝔽|ij=|ji`$. A Werner state is entangled iff $`\alpha (1,d]`$ and 1-distillable iff $`\alpha (d/2,d]`$. Moreover, it was shown in Iso  that any NPPT state can be mapped onto an entangled Werner state by means of LOCC operations. Therefore we can in the following restrict our discussion to the activation of Werner states keeping in mind that the obtained results hold for any NPPT state.

Isotropic states Iso ; Sym  commuting with all unitaries of the form $`U\overline{U}`$ (where $`\overline{U}`$ is the complex conjugate of $`U`$) are combinations of the maximally mixed state $`\mathrm{𝟏}/d^2`$ and the projector $`=|\mathrm{\Omega }\mathrm{\Omega }|`$ onto the maximally entangled state $`|\mathrm{\Omega }=1/\sqrt{d}_{i=1}^d|ii`$:

$$\omega (f)=f+\frac{1f}{d^21}\left(\mathrm{𝟏}\right),f[0,1].$$
(2)
An Isotropic state $`\omega `$ is known to be 1-distillable iff the maximally entangled fraction $`f=\mathrm{\Omega }|\omega |\mathrm{\Omega }>1/d`$, which is a sufficient condition for any other state as well Iso . Hence, an activation protocol succeeded if this condition is fulfilled by the final state.

The symmetric states playing the central role in the present paper act on a larger Hilbert space $`=_{A_1}_{A_2}_{B_1}_{B_2}`$ of total dimension $`d^4`$, where $`A`$ and $`B`$ again label the two parts of the system situated at different locations. The symmetry group under consideration is the group of all unitaries of the form $`W=(UV)_A(U\overline{V})_B`$. States $`\sigma `$ commuting with all these unitaries can most easily be expressed in terms of the minimal projectors $`\{P_i\}`$ spanning the commutant of the group commutant :

$`W:[\sigma ,W]=0`$ $``$ $`\sigma ={\displaystyle \underset{i=1}{\overset{4}{}}}\lambda _iP_i/\mathrm{tr}\left[P_i\right],`$ (3)
$`\text{with}P_{\genfrac{}{}{0pt}{}{1}{2}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\mathrm{𝟏}𝔽\right)_1_2,`$ (4)
$`P_{\genfrac{}{}{0pt}{}{3}{4}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\mathrm{𝟏}𝔽\right)_1\left(\mathrm{𝟏}\right)_2.`$ (5)
Note that, as labeled by the indices, the tensor products correspond to a split $`1|2`$ (and not $`A|B`$). Positivity and normalization of $`\sigma `$ requires $`\lambda _i0`$ and $`_{i=1}^4\lambda _i=1`$ such that any state of the considered symmetry can be characterized by a vector $`\stackrel{}{\lambda }^3`$ lying in a tetrahedron, which is given by these constraints. We note further that the set of symmetric states in (3) is abelian, i.e., all symmetric states commute with each other.

The activation protocol we use follows closely an idea of Ref. CDK . Initially the two parties $`A`$ and $`B`$ are supposed to share a Werner state $`\rho (\alpha )`$ acting on $`_0=_{A_0}_{B_0}`$ with $`dim_0=d^2`$ and a symmetric state $`\sigma `$ on $`_1_2`$ given by Eq.(3). After a local filtering operation is applied by projecting onto maximally entangled states $`_{A_{0,1}}`$ and $`_{B_{0,1}}`$ (acting on $`_{A_0}_{A_1}`$ respectively $`_{B_0}_{B_1}`$) the maximally entangled fraction of the resulting state on system 2 is given by

$$f(\rho (\alpha );\sigma ):=\frac{\mathrm{tr}\left[(\rho (\alpha )\sigma )(_{A_{0,1}}_{B_{0,1}}_2)\right]}{\mathrm{tr}\left[(\rho (\alpha )\sigma )(_{A_{0,1}}_{B_{0,1}}\mathrm{𝟏}_2)\right]}.$$
(6)
We know that $`\sigma `$ activates $`\rho (\alpha )`$ if $`f(\rho (\alpha );\sigma )>1/d`$. Since the output state of the protocol is itself isotropic, this condition is also necessary for the activation.

Of course, we are only interested in cases where $`\alpha (1,d/2]`$, i.e., $`\rho (\alpha )`$ is entangled but not 1-distillable, and $`\sigma `$ is in turn a PPT bound entangled state. The latter requires the classification of the symmetric states in (3), which is the content of the next section.

## III Classification and Activation

The following discussion will mainly take place in the three dimensional space given by the expansion coefficients $`\stackrel{}{\lambda }=(\lambda _1,\lambda _2,\lambda _3)`$ from Eq. (3). A vector $`\stackrel{}{\lambda }`$ corresponds to a (positive and normalized) symmetric state $`\sigma (\stackrel{}{\lambda })`$ iff $`\stackrel{}{\lambda }𝒮=\{\stackrel{}{v}^3|v_i0,_iv_i1\}`$, where the state space $`𝒮`$ is a tetrahedron.

The set $`𝒫`$ corresponding to normalized operators with a positive partial transpose can easily be obtained by observing that the symmetry group of the partially transposed operators $`\sigma ^{T_A}`$ in (3) is equal to the group of unitaries $`W`$ when interchanging the systems $`12`$. The respective minimal projectors $`\{Q_i\}`$ can therefore be obtained from the projectors $`\{P_i\}`$ simply by relabeling the systems $`12`$, and the $`k`$-th coordinate of an extreme point $`\stackrel{}{p}^{(i)}`$ of $`𝒫`$ is thus given by

$$p_k^{(i)}=\mathrm{tr}\left[Q_i^{T_A}P_k\right]/\mathrm{tr}\left[Q_i\right].$$
(7)
Hence, $`𝒫`$ is again a tetrahedron and Eq. (7) leads to the extreme points:

$$\begin{array}{ccc}\stackrel{}{p}^{(1)}=(\frac{d1}{2d},\frac{d1}{2d},\frac{1d^2}{2d})\hfill & ,& \stackrel{}{p}^{(3)}=(\frac{1}{2d},\frac{1}{2d},\frac{d+1}{2d}),\hfill \\ \stackrel{}{p}^{(2)}=(\frac{1d}{2d},\frac{1+d}{2d},\frac{(d1)^2}{2d})\hfill & ,& \stackrel{}{p}^{(4)}=(\frac{1}{2d},\frac{1}{2d},\frac{d1}{2d}).\hfill \end{array}$$

Straight forward linear algebra now allows us to compute the extreme points $`\{\stackrel{}{\tau }^{(i)}\}`$ of the intersection $`𝒮𝒫`$ corresponding to the set of symmetric PPT states, which is shown in Fig.1:

$$\begin{array}{ccc}\stackrel{}{\tau }^{(1)}=(0,0,\frac{1}{2})\hfill & ,& \stackrel{}{\tau }^{(2)}=(0,0,0),\hfill \\ \stackrel{}{\tau }^{(3)}=(\frac{1}{2d},\frac{1}{2d},\frac{d1}{2d})\hfill & ,& \stackrel{}{\tau }^{(4)}=(0,\frac{1}{d},0),\hfill \\ \stackrel{}{\tau }^{(5)}=(\frac{1}{d+2},0,0).\hfill & & \end{array}$$

The PPT state space given by the convex hull of these points can be divided into three parts: separable states and activating respectively not activating bound entangled states.

### III.1 Separable states

A vector $`\stackrel{}{\lambda }𝒮𝒫`$ corresponds to a separable symmetric state iff we can find any (not necessarily symmetric) separable state $`\rho _{sep}`$ such that $`\lambda _i=\mathrm{tr}\left[P_i\rho _{sep}\right]`$ (cf.Sym ).

A special case of a symmetric separable state is of course a tensor product of a Werner and an Isotropic state $`\sigma =\rho (\alpha )_1\omega (f)_2`$, where both are separable. In fact, if we choose these two states lying on the separable boundaries, i.e., $`\alpha \{d,1\}`$ and $`f\{0,1/d\}`$, we retrieve the extreme points $`\stackrel{}{\tau }^{(1)},\mathrm{},\stackrel{}{\tau }^{(4)}`$.

Another point $`\stackrel{}{\tau }^{(0)}=(\frac{1}{2d},0,\frac{d2}{4d})`$, that will also turn out to be an extreme point of the set of separable symmetric states, is obtained from the product state:

$$\rho _{sep}=|\mathrm{\Phi }^+\mathrm{\Phi }^+|_A|\mathrm{\Psi }^{}\mathrm{\Psi }^{}|_B,$$
(8)
where $`|\mathrm{\Phi }^+=\frac{1}{\sqrt{2}}\left(|11+|22\right)`$ and $`|\mathrm{\Psi }^{}=\frac{1}{\sqrt{2}}\left(|12|21\right)`$ are two-dimensional maximally entangled states.

The convex hull of the points $`\stackrel{}{\tau }^{(0)},\mathrm{},\stackrel{}{\tau }^{(4)}`$ (see Fig.1) already covers the entire separable region as we are going to show in the following that the complement of this polytope within $`𝒮𝒫`$ corresponds to bound entangled states.

### III.2 Bound entangled and activating

The equation $`f(\rho (\alpha );\sigma (\stackrel{}{\lambda }))=\frac{1}{d}`$ written out as $`_ic_i(\alpha )\lambda _i=0`$, with

$$c_i(\alpha )=\frac{\mathrm{tr}\left[(\rho (\alpha )P_i)(_{A_{0,1}}_{B_{0,1}}(d\mathrm{𝟏})_2)\right]}{\mathrm{tr}\left[P_i\right]}$$

is linear in $`\lambda _i`$ and thus defines a plane separating symmetric states activating $`\rho (\alpha )`$ from states apparently not activating it. The task is now to construct this separating plane depending on the parameter $`\alpha `$.

As we have already used above, the points $`\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)}`$ correspond to product states of the form $`\sigma =\rho (\alpha )_1\omega (\frac{1}{d})_2`$ for which $`f(\rho (\alpha );\sigma )=\frac{1}{d}`$ obviously holds for any Werner state $`\rho (\alpha )`$. Thus $`\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)}`$ are two fix points of the separating plane with respect to a variation of $`\alpha `$ and we need to know only one more point. For this purpose we consider the line $`\stackrel{}{\lambda }(t)=t\stackrel{}{\tau }^{(5)}+(1t)\stackrel{}{\tau }^{(1)}`$. Solving the equation $`f(\rho (\alpha );\sigma \left(\stackrel{}{\lambda }(t)\right))=\frac{1}{d}`$ yields the required third point with

$$t=\frac{2+d}{2\alpha +d}.$$
(9)
This is obviously a strict monotone function in $`\alpha `$ and it leads to the following properties of the separating plane:

1. For $`\alpha =\frac{d}{2}`$ corresponding to the boundary Werner state which is not 1-distillable, Eq. (9) leads to $`\stackrel{}{\lambda }(t)=\stackrel{}{\tau }^{(0)}`$, showing that $`\stackrel{}{\tau }^{(0)}`$ indeed lies on the boundary of the set of separable states. That is, any PPT state in front of the plane $`\{\stackrel{}{\tau }^{(0)},\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)}\}`$ must be bound entangled since it activates at least $`\rho (\alpha =\frac{d}{2})`$.
2. In the limit $`\alpha 1`$, i.e., $`\rho (\alpha )`$ becoming less and less entangled, $`\stackrel{}{\lambda }(t)`$ approaches $`\stackrel{}{\tau }^{(5)}`$. However, for any $`\alpha =1+\epsilon ,\epsilon >0`$ the polytope $`\{\stackrel{}{\tau }^{(0)},\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)},\stackrel{}{\lambda }(t)\}`$ has a nonempty interior corresponding to PPT bound entangled states capable of activating any $`\rho (\alpha )`$ with $`\alpha >1+ϵ`$.
3. Except for the line $`\{\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)}\}`$ all PPT states on the plane $`\{\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)},\stackrel{}{\tau }^{(5)}\}`$ lie on the activating side of the separating plane for any $`\alpha >1`$. The corresponding symmetric states can thus be considered to be universal activators in the sense that they activate any entangled Werner state and therefore any NPPT state.

The set of bound entangled universal activators contains states arbitrary close to the line $`\{\stackrel{}{\tau }^{(3)},\stackrel{}{\tau }^{(4)}\}`$ which in turn corresponds to separable states. By the continuity properties of the entanglement measures Entanglement of Formation EoF  and Relative Entropy of Entanglement RelEnt  this geometric vicinity, however, translates directly to the proposition that these states contain an arbitrary small amount of entanglement.

### III.3 Bound entangled and not activating

In order to complete the classification of the symmetric states introduced in Eq. (3) we have still to determine the entanglement properties of the states corresponding to the tetrahedron $`\{\stackrel{}{\tau }^{(0)},\stackrel{}{\tau }^{(2)},\stackrel{}{\tau }^{(4)},\stackrel{}{\tau }^{(5)}\}`$. The plane separating this set from the separable states derived above is characterized by a linear operator $`W`$ via $`\mathrm{tr}\left[W\sigma (\stackrel{}{\lambda })\right]=0`$, where

$$W=\left(\mathrm{𝟏}𝔽\right)_1\left(\mathrm{𝟏}\frac{d}{2}\right)_2.$$
(10)
However, this operator is an entanglement witness (cf.optwit ), meaning that $`\mathrm{tr}\left[W\rho \right]0`$ holds for any separable state $`\rho `$. In order to see this property we have just to utilize the results from IHa ; IBM , where it was shown that $`\mathrm{𝟏}\frac{d}{2}`$ has a positive expectation value on any pure state of Schmidt rank two. Since the antisymmetric projector $`\frac{1}{2}\left(\mathrm{𝟏}𝔽\right)`$ is a sum of Schmidt rank two states, $`\varphi _A\psi _B|W|\varphi _A\psi _B`$ is a sum of such positive expectations for any pure product state $`|\varphi _A\psi _B`$. Hence, $`\mathrm{tr}\left[W\rho \right]0`$ is indeed fulfilled by any separable state implying that the tetrahedron under discussion corresponds to bound entangled symmetric states which are however not activating (with respect to the considered protocol).

## IV Conclusion

We investigated the entanglement properties of a new abelian set of symmetric states with regard to the activation of NPPT distillation. The set contains PPT bound entangled states of full rank providing a universal resource for the activation of any NPPT state (in any finite dimension). Some of these universal activators lie arbitraryly close to separable states such that the activation process turns out to require only an infinitesimal amount of entanglement. This indicates that the difference between 1-distillable and $`n`$-undistillable or even bound entangled NPPT states is very subtle. Moreover, the above results show that even weakly entangled bound entangled PPT states can be useful for some quantum information processing purposes BEtele .

As the problem discussed in this paper is primarily a feasibility problem, we have at this stage not asked for the obtained rates. In fact, one could for instance easily improve the probability of success for the used activation protocol by a factor $`d^2`$ by measuring in a basis of maximally entangled states and retaining the state whenever the measurement outcomes coincide. An interesting question going one step further and requiring knowledge about rates is, whether there is the possibility of self-activation after some initial activation with a limited resource took place. In other words, is it possible to yield asymptotically more entanglement from distillation than is needed for the preparation of the activator states?

## Acknowledgement

The authors would like to thank R.F. Werner and T. Eggeling for stimulating discussions and in particular M. Lewenstein for bringing the witness property of $`W`$ to their attention. Funding by the European Union project EQUIP (contract IST-1999-11053) and financial support from the DFG (Bonn) is gratefully acknowledged.