text:
# The 3He(e, e′d)p Reaction in (𝑞,𝜔)-constant Kinematics

## 1 Introduction

Many nuclear properties can be described successfully within a mean-field approach. However, phenomena like the depletion of spectroscopic strength and the occurrence of bumps at missing energies characteristic of two-nucleon emission in $`(e,e^{}p)`$ reactions, indicate that correlations between nucleons, i.e., the motion inside a nucleus of two nucleons as a pair with a certain relative motion, also play an essential role. The ($`e`$,$`e^{}d`$) reaction has proven to be a sensitive tool for the investigation of proton-neutron ($`pn`$) correlations in nuclei . Assuming direct knock-out of a $`pn`$ pair, the cross section for the ($`e`$,$`e^{}d`$) reaction can approximately be written  as $`d^6\sigma /dE_e^{}d\mathrm{\Omega }_e^{}dE_dd\mathrm{\Omega }_d=K\sigma _{e,pn}(q)S_{pn}(E_m,p_m,p_d)`$. Here, the cross section $`\sigma _{e,pn}`$ for scattering of an electron from a $`pn`$ pair leading to a deuteron in the final state, which depends on the momentum transfer $`q`$, reflects the relative proton-neutron motion, i.e., the relative $`pn`$ wave function. The distorted spectral function $`S_{pn}`$, which depends on the missing energy $`E_m`$, the missing momentum $`p_m`$ and the momentum of the final deuteron $`p_d`$, contains the information about the centre-of-mass (c.o.m.) motion of the $`pn`$ pair within the nucleus, modified by final-state interactions (FSI). Since the deuteron has isospin zero and the initial $`pn`$ pair can be in a $`T=0`$ or $`T=1`$ state, in general both $`\mathrm{\Delta }T=0`$ and $`\mathrm{\Delta }T=1`$ transitions are possible. The former resembles quasi-elastic knockout of a deuteron, whereas the latter is similar to the (inverse of) deuteron disintegration. Both types of transitions have been studied previously (, ).

In this paper we present data for the $`^3`$He$`(e,e^{}d)`$ reaction, taken in so-called ($`q,\omega `$)-constant kinematics, in which the energy and momentum transfer to the nucleus are held constant, while the angle of the outgoing deuteron is varied. The data is compared to the results of three-body calculations  and to those of a covariant and gauge-invariant diagrammatic approach . For the three-nucleon system Faddeev calculations are nowadays available both for the ground state and for the continuum. The first available $`^3`$He$`(e,e^{}d)`$ data  has been compared to the results of such calculations . The emphasis in these calculations was on the correct treatment of FSI, which turned out to be crucial, especially for kinematics dominated by deuteron knock-out processes. The same data was also employed in Refs. , with the emphasis on the consistency between the dynamics and the one- and many-body currents. The agreement with the data, especially for parallel kinematics, was in both cases fairly good. Other data  was of limited kinematic coverage and no comparison to Faddeev calculations was made. More recently, we published data for the $`^3`$He$`(e,e^{}d)`$ reaction taken in parallel kinematics . Theoretical calculations gave a fair description of the data, but there were some discrepancies, especially for the $`q`$ dependence, and for the transverse structure function $`W_T`$.

Data in ($`q,\omega `$)-constant kinematics extend the range in which the theoretical predictions can be tested, and feature some new aspects. First of all one can investigate the longitudinal-transverse (LT) interference structure function, which is more sensitive to details of the interaction currents. Furthermore, one can investigate the dependence of the longitudinal (L) and transverse (T) structure functions on the kinematics. Finally, in selected parts of this kinematics one gets non-negligible contributions from processes in which the photon hits the non-detected proton. Interference with the process where the photon interacts with (the nucleons of) the deuteron gives an extra sensitivity to a correct description of the reaction.

Since the publication of the data obtained in parallel kinematics there have been several theoretical developments (see section 3). The Faddeev calculations by Golak et al.  were extended by using the AV18 NN potential in addition to the Bonn-B potential, and by including meson-exchange currents (MEC). Furthermore, the number of partial waves was increased to ensure full convergence. New calculations in the diagrammatic approach by Nagorny  employed also the AV18 potential instead of the Reid Soft Core potential, and now include all angular momentum states in the $`^3`$He ground state.

Since the data to be discussed comprise both the $`q`$ and $`p_m`$ dependence of the reaction, different aspects of the $`pn`$ motion in the nucleus 3He can be studied. Furthermore, the cross sections were measured at two beam energies, so that a separation of the longitudinal ($`W_L`$) plus transverse-transverse ($`W_{TT}`$), transverse ($`W_T`$), and longitudinal-transverse interference ($`W_{LT}`$) structure functions could be performed.
For the definition of the structure functions we follow Raskin and Donnelly , who write the differential cross section for the unpolarized ($`e`$,$`e^{}d`$) reaction as:

$$\frac{d^5\sigma }{dE_e^{}d\mathrm{\Omega }_e^{}d\mathrm{\Omega }_d}=C\left(v_LW_L+v_TW_T+v_{LT}W_{LT}\mathrm{cos}\varphi +v_{TT}W_{TT}\mathrm{cos}2\varphi \right),$$
(1)
where $`\varphi `$ is the angle between the electron scattering plane and the plane defined by the momentum transfer $`\stackrel{}{q}`$ and the momentum of the outgoing deuteron $`\stackrel{}{p_d}`$. The factor $`C`$ contains the Mott cross section and kinematical factors. Its precise form, as well as that of the kinematical factors $`v_i`$, are given in Ref. . The structure functions $`W_{LT}`$ and $`W_{TT}`$ are zero in parallel kinematics.

The separate structure functions, which result from different combinations of the components of the nuclear current, have a different sensitivity to the various aspects of the reaction. For instance, the coupling of the virtual photon to a $`T=1`$ $`pn`$ pair, transforming it into a deuteron, involves a spin-flip and thus is purely transverse, whereas the coupling to an initial $`T=0`$ pair, which resembles elastic $`ed`$ scattering, is dominantly longitudinal at our values of $`q`$ (see also ). Furthermore, one expects that MEC will mainly contribute to $`W_T`$. On the other hand, according to the three-body calculations the effects of FSI on $`W_L`$ are quite different from the ones on $`W_T`$.

## 2 Experiment

The $`^3`$He$`(e,e^{}d)p`$ reaction has been measured “in-plane”, using ($`q,\omega `$)-constant kinematics, for missing momenta up to 210 MeV/$`c`$, and for values of the transferred momentum $`q`$ of 412, 504 and 604 MeV/$`c`$. The (central) kinematics is given in Table 1.

In later figures a negative sign of $`p_m`$ denotes kinematics in which $`\theta _d<\theta _q`$ (corresponding to $`\varphi =0`$ in Eq. 1), and a positive $`p_m`$ kinematics with $`\theta _d>\theta _q`$ (corresponding to $`\varphi =\pi `$ in Eq. 1).
Since the experiment has been described in detail in Refs. , only the main points are mentioned here.

The experiment was performed with the extracted electron beam from the pulse-stretcher ring AmPS  at NIKHEF. The beam energies were 370 and 576 MeV. The extracted electron current was about 6 $`\mu `$A and had a duty factor of about 70%. The scattered electrons were detected with the QDD spectrometer and the knocked-out deuterons with the QDQ spectrometer .

A cryogenic gas target  operating at 20 K and 1.5 Mpa was used, which was filled with a mixture of 3He and 4He gases. In this way data was collected simultaneously for the three reactions: $`^3`$He$`(e,e^{}d)p`$, $`^4`$He$`(e,e^{}d)d`$ and $`^4`$He$`(e,e^{}d)pn`$. Results on the latter reaction channel have been published separately .

The resolution in missing energy varied between 0.3 and 2.0 MeV, depending on the kinematics, which was sufficient to separate the different reaction channels. The absolute 3He and 4He target thicknesses were determined by comparison of the measured elastic scattering cross sections to calculated ones . During the ($`e`$,$`e^{}d`$) measurements the total target thickness was monitored through the singles rate of either one of the spectrometers. Checks with elastic scattering before and after each set of ($`e`$,$`e^{}d`$) measurements were consistent to within 2%.

Both the 2H($`e`$,$`e^{}d`$) and the 1H($`e`$,$`e^{}p`$) reaction were used in order to check the reconstruction of the electron and hadron momentum vectors in the QDD and QDQ spectrometers and to determine the coincidence detection efficiency.

The data analysis included the following steps. First, the deuterons were separated from protons and tritons by using the pulse height from the scintillators in the QDQ spectrometer. Next, the particle vectors at the target were reconstructed. Since an extended target was used, this reconstruction and the acceptances of both spectrometers depend on the position of the interaction point along the beam. By using energy and momentum conservation the values of the missing energy

$$E_m=\omega T_pT_d,$$
(2)
where $`T_x`$ is the kinetic energy of particle $`x`$, and the missing momentum

$$p_m=\stackrel{}{q}\stackrel{}{p_d}$$
(3)
were calculated from the particle vectors. Next, the accidental coincidences were subtracted. As a result of the high duty-factor of the extracted beam, the real-to-accidental ratio was high, ranging from 15 to 1640 depending on the kinematics. Next, the data was normalised to the target thickness, the integrated charge and the experimental detection volume, and corrected for detection inefficiencies. The detection volume was obtained from a Monte-Carlo simulation, which uses the measured optical properties of the spectrometers , including their vertex-position dependent angular acceptances. Finally, the data was radiatively unfolded. The systematic uncertainty amounts to about 3% for the cross sections and 4-5% for the structure functions .

The cross section measured in one kinematical setting covers an appreciable range in $`p_m`$, as a consequence of the angular and momentum acceptances of the spectrometers. However, over this range the values of $`q`$ and the kinematical factors in the cross-section expression of Eq. (1) vary slightly around the central values. A Monte-Carlo simulation was performed to determine the average values of the kinematical factors and of $`q`$ for the different $`p_m`$ bins within one measurement. Then, the cross sections for the different $`p_m`$ bins were recalculated to a common $`q`$ value by using the $`q`$ dependence of the cross section as measured in this experiment. The correction to the cross section amounted to typically a few percent.

## 3 Theoretical calculations

The experimental data is compared to the calculations of Van Meijgaard and Tjon , Golak et al.  and Nagorny et al. . The results of Van Meijgaard and Tjon are based on solutions of the Faddeev equations for the three-body system, employing a central local $`NN`$-interaction, the spin-dependent Malfliet-Tjon I-III potential , in the Unitary Pole Expansion (UPE). Since this interaction contains only $`s`$-wave forces, the ground-state wave function of 3He includes only $`s`$-waves. Furthermore, a relativistic current operator is used.

The calculations of Golak et al. are also based on the solution of Faddeev-like equations, but employ the one-boson-exchange Bonn-B potential  and the AV18 interaction. . Thus state-of-the-art $`^3`$He and 3N continuum wave functions are employed. A non-relativistic single-nucleon current operator is used, consistent with the non-relativistic Schrödinger equation for the 3N states. (The lowest-order relativistic corrections to the single-nucleon density operator were studied and found to be on the $`5\%`$ level in the range of small $`p_m`$ values.) In addition $`\pi `$\- and $`\rho `$-like exchange currents  are included, which in case of the AV18 potential are consistent with the forces. The Riska prescription  was employed, which guarantees that the currents fulfill the continuity equation. In this approach the transverse currents $`j_x`$ and $`j_y`$ explicitly include MEC, while $`j_z`$ is related to the charge current $`j_0`$ through the continuity equation. MEC were not included in the charge density operator. This would be a relativistic effect, which was not considered. In case of the Bonn-B potential standard  $`\pi `$\- and $`\rho `$-like exchange currents are taken. The electromagnetic nucleon form factors used are the ones from Hoehler et al. .

Nagorny et al. use quite a different approach, which is based on including the electromagnetic field into the strongly-interacting system in a fully relativistic and gauge-invariant way . Two covariant sets of diagrams, including pole, “contact” and one-loop diagrams are used, which provide both nuclear-current conservation and inclusion of the dominant part of FSI and MEC effects in a form that is consistent with the nuclear dynamics. The diagrams were generated by “minimal insertion” of the electromagnetic field into all external/internal lines and also directly into the 3- and 4-point nuclear vertices, which produces various contact currents in accordance with Ward-Takahashi identities . The strong form factors in the covariant nuclear vertices 3He$`pd`$ and 3He$`ppn`$ are taken as the positive-energy states in the laboratory frame through the solutions of the Faddeev equations with the AV18 potential. All angular momentum states in the ground state of $`^3`$He are included. The electromagnetic form factors in the completely relativistic $`\gamma NN`$-, $`\gamma dd`$\- and $`\gamma ^3`$He3He-vertices are taken from standard parameterizations of experimental form factors.

## 4 Results

### 4.1 Direct proton knock-out

We first focus on the measurements at $`q=412`$ MeV/$`c`$ at the higher beam energy, where the largest range in $`p_m`$ has been probed. The measured $`^3`$He$`(e,e^{}d)`$ cross sections for this kinematics are shown in Fig. 1.

Above a certain value of $`p_m`$ the cross section flattens off, an effect that is not seen in parallel kinematics . The origin of this effect can easily be understood. If one neglects FSI the virtual photon can either couple to the proton that is knocked out, while the detected deuteron acts as a spectator (direct proton knock-out: DPKO), or the virtual photon couples to either one of the nucleons that constitute the final deuteron, which is knocked out, while the proton is a spectator (direct deuteron knock-out: DDKO). Inclusion of FSI modifies this picture quantitatively, but not qualitatively.

Now one can look into the kinematics. In the case of DDKO, $`p_m`$ and (the constant value of) $`q`$ combine to the detected deuteron momentum $`p_d`$, which decreases from about 400 MeV/$`c`$ at $`p_m=0`$ MeV/$`c`$ to about 340 MeV/$`c`$ at $`p_m=200`$ MeV/$`c`$. However, in case of DPKO the effective $`p_m`$ is equal to $`p_d`$. Since for both processes the cross section decreases with the effective value of $`p_m`$, the cross section is dominated at low $`p_m`$ by DDKO, but for increasing $`p_m`$ values the contribution of DPKO increases.

Although for a quantitative comparison with the data the effects of distortions have to be included, the PWIA and BORN calculations of Van Meijgaard and Tjon show this nicely. The PWIA calculation contains only the DPKO term, while the BORN calculation includes both the DPKO and the DDKO contribution. Whereas for $`p_m=0`$ the PWIA calculation is orders of magnitude below the data, it increases rapidly with $`p_m`$, ending up only one order of magnitude below the data at the largest values of $`p_m`$ measured. The BORN calculation is largest around $`p_m=0`$ and then decreases, but flattens off at large $`p_m`$ values due to the increasing contribution of the DPKO part.

### 4.2 Cross sections

Before discussing all data and the comparison with the calculations, it is of interest to look into the influence of final-state interactions (FSI). The ratio of the FULL cross sections, i.e., including FSI, over the ones calculated in PWIAS, which means a fully antisymmetrized plane-wave impulse approximation  (equivalent to the BORN cross section by Van Meijgaard and Tjon), is shown in Fig. 2. One observes that FSI effects around $`p_m=0`$ decrease the cross section considerably at $`q=412`$ MeV/$`c`$, but less so at larger $`q`$. This has already been noticed before . Here, one should keep in mind that the $`pd`$ centre-of-mass energy is only 14.7 MeV at the lowest $`q`$-value, rising to 31.1 MeV at the highest $`q`$-value. In the latter case the reduction of the cross section at low values of $`p_m`$ is already rather small, about 10-15%. Beyond a $`|p_m|`$-value of about 150 MeV/$`c`$ FSI increases the cross section considerably. (Overall the effects are slightly smaller at the lower electron energy. We will come back to this later when discussing the separated structure functions).

The measured ($`e`$,$`e^{}d`$) cross sections for the various kinematics are shown in Fig. 3. Since the calculations of Golak et al. are the most detailed, we will first compare those to the data. One can make three observations. The first one is that the calculations overestimate the data. This holds already for the results with the Bonn-B potential, and even more for the AV18 ones. The inclusion of MEC increases the discrepancy. Secondly, the (relative) difference between data and calculations is slightly smaller at the lower beam energy (we will discuss this further in the next subsection) and increases with $`q`$. Finally, the $`p_m`$ dependence of the data is fairly well described, except at $`q=412`$ MeV/$`c`$, where there is a discrepancy at negative $`p_m`$, in the region where the cross section flattens off due to the influence of DPKO. Obviously the interference between the DPKO and DDKO processes is not well described.

In Fig. 4 the three theoretical approaches are compared with each other and with the data. For the calculations by Golak et al. the AV18 results are chosen, since these forces and the accompanying MEC are considered to be the most realistic ones. In addition to what has been said already about the various calculations some more points are worth mentioning.

It is seen, especially at small $`p_m`$-values, that with increasing values of $`q`$ the calculations by Van Meijgaard and Tjon predict increasingly smaller values for the cross section compared to those calculated by Golak et al.. Since in first approximation the $`q`$-dependence of the cross section probes the wave function of the $`pn`$-pair, part of this may reflect the increasing importance of the $`D`$-state of the $`pn`$-pair (compare elastic scattering from the deuteron at these $`q`$-values ). As mentioned, the calculations by Van Meijgaard and Tjon do not contain $`d`$-waves. Golak et al. have also performed calculations with the MT I-III potential, using their non-relativistic single-nucleon current. The results agree quite well with the ones by Van Meijgaard and Tjon.

The calculations by Nagorny give a fairly good description of the data. At $`q`$ = 412 MeV/$`c`$ earlier calculations  were significantly lower than the data for $`p_m>`$100 MeV/$`c`$ (see the dash-dotted line in Fig. 4). In contrast to those the present calculations include also the $`l=2,L=2`$ configuration in the $`^3`$He ground state. Although this component is small, it has a relatively large influence at larger values of $`p_m`$ due to its overlap with the deuteron $`l=2`$ (D) state. Altogether, these calculations give a fairly good description of both the $`p_m`$\- and the $`q`$-dependence of the data.

### 4.3 Structure functions

As mentioned in section 2, the experimental data was obtained for values of $`\varphi `$ close to 0 and $`\pi `$, and for two values of the incoming electron energy, so the structure functions $`W_L+\frac{v_{TT}}{v_L}W_{TT}`$, $`W_T`$ and $`W_{LT}`$ (see Eq. 1) could be separated.

Before discussing the data it is instructive to look at the theoretical predictions. In Fig. 5 the separate structure functions at $`q=412`$ MeV/$`c`$ calculated by Golak with the AV18 interaction, with and without inclusion of FSI and MEC, are shown.

It is clear that $`W_L`$ and $`W_T`$ are the dominant structure functions, and that $`W_T`$ is not small compared to $`W_L`$. This is due to the contribution of the (fully transverse) $`\mathrm{\Delta }T=1`$ transition (the $`\mathrm{\Delta }T=0`$ transition is dominantly longitudinal), see Refs. . The $`W_{LT}`$ interference structure function is rather small, while $`W_{TT}`$ is negligible except for values of $`p_m`$ above about 150 MeV/$`c`$.
Another interesting observation is that the influence of FSI is large in $`W_L`$ (and in $`W_{LT}`$ and $`W_{TT}`$), but very small in $`W_T`$, except at large values of $`p_m`$. In this context one should realize that a PWIAS calculation for $`W_L`$ violates current conservation. On the other hand MEC, as included in the present approach, only influence the ’transverse current’ dependent $`W`$’s, considerably increasing the values of $`W_T`$ and $`W_{TT}`$.

$`W_{LT}`$

The $`W_{LT}`$ structure function can be determined by just comparing the data in kinematics with $`\theta _d`$ smaller and larger than $`\theta _q`$ (see section 2). The separated $`W_{LT}`$ structure functions for the three values of $`q`$ are shown together with the calculations in Fig. 6, both for $`E_e`$= 576 MeV and $`E_e`$= 370 MeV. The values obtained at the two different energies agree within the experimental uncertainties. As is already clear on the cross section level, where the difference between the data for $`\theta _d`$ smaller and larger than $`\theta _q`$ is not very large, the measured $`W_{LT}`$ is small, so the error bars are relatively large. The theoretical curves for $`W_{LT}`$ were obtained by separating the calculated cross sections in the same way as the measured cross sections. (Naturally they agreed with the ones calculated directly, using the appropriate components of the hadronic current operator). For all calculations $`W_{LT}`$ is independent of the incoming energy, as expected, and has a maximum near $`p_m`$= 50 MeV/$`c`$ for all $`q`$ values, while it decreases for higher $`p_m`$ values. In the Faddeev calculations FSI has two effects on $`W_{LT}`$: both the magnitude of the maximum decreases and the decrease as a function of $`p_m`$ is stronger. This even introduces a zero crossing at $`q`$= 504 MeV/$`c`$ and a double zero crossing at $`q`$= 412 MeV/$`c`$. The latter presumably is related to the increasing influence of DPKO, which is also visible in this region in the cross section.
In the $`p_m`$ range between 0 and about 70 MeV/$`c`$ the data and the various calculations agree with each other within the experimental accuracy. However, at larger values of $`p_m`$ the calculations tend to underestimate the data.

$`W_L`$ and $`W_T`$

Data at $`q`$= 412 MeV/$`c`$ and $`q`$= 504 MeV/$`c`$ were taken at both energies. Therefore, for the first time the $`W_L(+\frac{v_{TT}}{v_L}W_{TT})`$ and $`W_T`$ structure functions could be separated in ($`q,\omega `$) kinematics, for those $`p_m`$ values where all four cross sections were available. Since all calculations predict $`W_{TT}`$ to be very small compared to $`W_L`$, except at values of $`p_m`$ above about 150 MeV/$`c`$, this is effectively an L-T separation. The results are presented in Fig. 7. It turns out that the values found for $`W_L`$ and $`W_T`$ in ($`q,\omega `$) kinematics are very close to the values found in parallel kinematics . This, together with the fact that $`W_{LT}`$ and $`W_{TT}`$ are small, suggests that a factorization of the cross section into an ’$`epn`$’ cross section and a (distorted) spectral function (see also the introduction and Ref. ) is a reasonable approximation.

The calculations by Golak et al. with the AV18 potential and including MEC are in global agreement with $`W_T`$, but overestimate $`W_L`$, especially at $`q=504`$ MeV/$`c`$. In this context we recall that at $`q=604`$ MeV/$`c`$, where no separation of the structure functions could be performed, the unseparated cross section was even more overestimated. Since the $`q`$-dependence in first instance probes the relative motion of the $`pn`$ pair in $`^3`$He, this suggests that the short-distance behaviour of this motion is not well described. (As the used current operator is a non-relativistic one, it was verified that elastic scattering from the free deuteron was well described).

Whereas the calculations by Nagorny gave a fairly good description of the measured cross sections, they clearly underestimate the longitudinal structure function $`W_L`$ at $`q=412`$ MeV/$`c`$, and overestimate $`W_T`$. The calculations show a minimum for $`W_L`$, which is not present in a PWIAS calculation. This has been observed before by Nagorny et al. . The cause is a large, destructively interfering, contribution of the 3He-pole diagram (S-term) in ($`q,\omega `$) kinematics. This term, which is part of the FSI, is constant with $`p_m`$, since it only depends on $`\omega `$. Another part of the FSI comes from the “contact current”, but this part is still small at $`q=412`$ MeV/$`c`$. Clearly, this indicates an inadequacy of the treatment of FSI. Presumably, inclusion of only the pole part of $`pdpd`$ rescattering is not sufficient, and the regular part  has to be included as well.

Although the pioneering calculations by Van Meijgaard and Tjon give a reasonable description of the data, this observation does not allow strong conclusions, since the used interaction is rather simple, no MEC are included, and only $`s`$-waves are taken into account.

## 5 Summary and conclusions

Cross sections for the $`^3`$He$`(e,e^{}d)p`$ reaction have been measured in ($`q,\omega `$)-constant kinematics for a range in missing momentum at two beam energies and for three values of the three-momentum transfer. Thus the separate L(+TT), T and LT structure functions could be determined at $`q`$ = 412 and 504 MeV/$`c`$, which provide for a sensitive test of theoretical calculations for this reaction. The LT structure function is found to be rather small, but the transverse structure function is of comparable magnitude as the longitudinal one. This points to the importance of the $`\mathrm{\Delta }T=1`$ transition, in which a $`T=1`$ $`pn`$ pair in $`^3`$He is transformed into a $`T=0`$ deuteron in the final state.

The data has been compared to the results of two types of three-body Faddeev calculations, one with a simple interaction and containing only $`s`$-waves, the other one employing the Bonn-B and the AV18 potentials, and including meson-exchange currents (MEC). Also calculations based on a covariant diagrammatic approach, including tree and one-loop diagrams, and using Ward-Takahashi identities to ensure gauge invariance were used.

All calculations give a fair description of the $`p_m`$\- and the $`q`$-dependences of the cross sections, which reflect the centre-of-mass and relative motion of the $`pn`$ pair inside $`^3`$He, respectively, and of their L/T character. This means that the essential ingredients of the calculations: the structure of $`^3`$He in its ground state, and final-state interactions (the continuum structure of the 3N system), are reasonably well understood. Final-state interactions have a large influence, especially at the lower value of $`q`$ (and accompanying low centre-of-mass energy of the final $`p`$-$`d`$ system). At low values of $`p_m`$ they reduce the longitudinal part of the cross section, whereas at large $`p_m`$-values they increase it. In this region also direct proton knock-out starts to noticeably influence the cross sections, leading to a rise of the cross section. This was experimentally observed at $`p_m`$-values above 150 MeV/$`c`$ at the lowest value of $`q`$. In contrast the value of the transverse structure function hardly changes when including FSI. Meson-exchange currents increase considerably the transverse structure function. All calculations predict relatively small values for the LT interference structure function, which is consistent with the data. Also $`W_{TT}`$ is calculated to be small, but this structure function cannot be separated with the present “in-plane” data.

However, upon closer inspection there are also significant differences between the data and the theoretical calculations. The diagrammatic approach gives a fair description of the cross section, but under(over)estimates the longitudinal (transverse) structure function. It also predicts a minimum in the longitudinal structure function, presumably due to the neglect of the regular part of $`pdpd`$ rescattering. The most striking result is that the supposedly best Faddeev calculations available at present, which employ the AV18 nucleon-nucleon interaction and include MEC, overestimate the measured cross sections, the more the larger the value of $`q`$. The major discrepancy is in the longitudinal cross section, the transverse one being reasonably well described. Possible explanations for this discrepancy could be: the use of the Sachs form factors $`G_E^p`$ and $`G_E^n`$ instead of the Pauli form factors $`F_1^p`$ and $`F_1^n`$ reduces the longitudinal response function $`R_L`$ by about 15-25%, depending on the kinematics. Furthermore, inclusion of a three-nucleon force will influence the binding energy and hence the $`dp`$ overlap. A PWIAS calculation including the Urbana IX 3N force  for the $`^3`$He bound state gives a reduction of the cross section of about 15%. Finally meson-exchange currents should also be included in the charge-density operator, since it is known that they influence e.g. the elastic charge form factor. However, such calculations are not yet available.

This work is part of the research programme of the “Stichting voor Fundamenteel Onderzoek der Materie (FOM)”, which is financially supported by the “Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)”. The Faddeev calculations by Golak et al. were supported by the Deutsche Forschungsgemeinschaft and the Polish Committee for Scientific Research. The calculations were performed on the T90 of the NIC in Jülich.


format: text

text:
# Stable spinning optical solitons in three dimensions

## Abstract

We introduce spatiotemporal spinning solitons (vortex tori) of the three-dimensional nonlinear Schrödinger equation with focusing cubic and defocusing quintic nonlinearities. The first ever found completely stable spatiotemporal vortex solitons are demonstrated. A general conclusion is that stable spinning solitons are possible as a result of competition between focusing and defocusing nonlinearities.

Optical solitons (spatial, temporal, or spatiotemporal) are self-trapped light beams or pulses that are supported by a balance between diffraction and/or dispersion and various nonlinearities. They are ubiquitous objects in optical media George . Spatiotemporal solitons (STS) KanRub , alias superspikes MBA ; Manassah  or light bullets Yaron , were found in many works KanRub  \- Agrawal . Although they cannot be stable in the uniform Kerr ($`\chi ^{(3)}`$) medium collapse , stability can be achieved in saturable MBA ; Edm , quadratically nonlinear ($`\chi ^{(2)}`$) KanRub ; quadr , and graded-index Kerr media Agrawal . While a fully localized STS in three dimensions (3D) has not yet been found in an experiment, 2D ones were observed in a bulk $`\chi ^{(2)}`$ medium Wise . The interplay of spatio-temporal coupling and nonlinearity may also play an important role in self-defocusing media ManassahDefocusing .

Spinning (vortex) solitons are also possible in optical media. Starting with the works vortex , both delocalized (“dark”) and localized (“bright”) optical vortices in 2D were investigated DiTrapp ; unstable1 ; experiment . In the 3D case they take the shape of a torus (“doughnut”) Anton ; new . However, the only previously known physical model which could support stable 3D vortex solitons is the Skyrme model reviews , which has recently found a new important application to Bose-Einstein condensates (BEC) SkyrmeBEC . Our objective in this paper is to identify fundamental models of the nonlinear-Schrödinger (NLS) type in 3D that give rise to stable spinning solitons, as NLS models are much simpler and closer to more experimental situations, having applications to optics, BEC, plasmas, etc. (see below).

For bright vortex solitons stability is a major issue as, unlike their zero-spin counterparts, the spinning solitons are prone to destabilization by azimuthal perturbations. In 2D models with $`\chi ^{(2)}`$ and saturable nonlinearities an azimuthal instability was revealed by simulations unstable1  and observed experimentally experiment . As a result, a soliton with spin $`1`$ splits into two or three fragments, each being a moving zero-spin soliton. Simulations of the 3D spinning STS in the $`\chi ^{(2)}`$ model also demonstrates its instability-induced splitting into separating zero-spin solitons new . Nevertheless, the $`\chi ^{(2)}`$ nonlinearity acting in combination with the self-defocusing Kerr ($`\chi ^{(3)}`$) nonlinearity, gives rise to the first examples of stable spinning (ring-shaped) 2D solitons with spin $`s=1`$ and $`2`$ Isaac . It should be stressed that all the 2D spinning solitons actually represent static spatial beams; on the contrary, 3D solitons are moving spatiotemporal ones, which are localized not only in the transverse plane, but also in the propagation coordinate, see below.

A model which may support stable spinning solitons in 3D is the one with a cubic-quintic (CQ) nonlinearity, which (in terms of optics) assumes a nonlinear correction to the medium’s refractive index in the form $`\delta n=n_2In_4I^2`$, $`I`$ being the light intensity. The CQ nonlinearity was derived, starting from the Maxwell-Bloch equations, for light propagation combining resonant interaction with two-level atoms and dipole interactions between the atoms dipole , or the Kerr nonlinearity of a waveguide John  (see also Ref. Manassah-twolevel ). A unifying feature of those media is competition between different nonlinearities, and stable solitons may exist in the range of intensities where the competition takes place (which may be controlled, for instance, through the density of two-level atoms in a $`\chi ^{(3)}`$ waveguide). In fact, the NLS equation of the CQ type is a generic model, which also applies to Langmuir waves in plasmasplasma  and BEC BEC  (although in the latter case, three-body interactions, which give rise to the quintic term, may induce losses through recombination of the colliding atoms, thus making the quintic coefficient complex).

In the first simulations of 2D solitons with spin $`1`$ in the CQ model, it was found that they propagated in a stable way, provided that their energy is not too small Q . A later analysis, based on the computation of linear-stability eigenvalues, demonstrated that some of the spinning 2D solitons considered in Ref. Q  are subject to a weak azimuthal instability. Nonetheless, in another part of their existence region, with very large energies, solitons with spin $`s=1`$ and $`s=2`$ were confirmed to be truly stable in the 2D CQ model IsaacPLA  (all the solitons with $`s3`$ are unstable).

It was recently shown by direct simulations of the CQ model new2  that 3D spinning solitons with moderate energies are unstable against azimuthal perturbations, while the ones with very large energies, i.e., broad “doughnuts” with a small hole in the center, were robust under propagation. However, a consistent stability analysis makes it necessary to compute eigenvalues of small perturbations. We will conclude that sufficiently broad STS with spin $`s=1`$ are stable, the stability region occupying $`20\%`$ of their existence region, while all the STS with $`s2`$ are unstable.

The evolution of the electromagnetic field envelope $`A`$ in the dispersive CQ medium is governed by the NLS equation,

$$2i\kappa _0A_z+_{}^2A+\kappa _0DA_{\tau \tau }+2\kappa _0^2(n_2/n_0)|A|^2A2\kappa _0^2(n_4/n_0)|A|^4A=0,\tau tz/V,$$
(1)
where $`z`$ and $`t`$ are the propagation coordinate and time, $`\kappa _0`$ and $`V`$ are the propagation constant and group velocity of the carrier wave, $`D>0`$ is the temporal dispersion, $`_{}^2`$ acting on the transverse coordinates $`x`$ and $`y`$. Equation (1) does not include higher-order effects, such as self-steepening, stimulated Raman scattering, non-paraxial diffraction, and third-order temporal dispersion (which, in another context, were taken into regard for spatiotemporal superspikes in Refs. Manassah ), as we anticipate that only broad solitons (with the temporal width $`1`$ ps), for which these effects are small, may be stable.

Defining rescaled variables $`u=\sqrt{n_4/n_2}A`$, $`T=n_2\sqrt{2\kappa _0/Dn_0n_4}\tau `$, $`Z=\left(\kappa _0n_2^2/n_0n_4\right)z`$, and $`(X,Y)=\kappa _0n_2\sqrt{2/n_0n_4}(x,y)`$, we transform Eq. (1) into a normalized form new2 ,

$$iu_Z+\left(u_{XX}+u_{YY}+u_{TT}\right)+|u|^2u|u|^4u=0.$$
(2)
STS solutions to Eq. (2) are sought as $`u=U(r,T)\mathrm{exp}(is\theta )\mathrm{exp}(i\kappa Z)`$, where $`r`$ and $`\theta `$ are the polar coordinates in the transverse plane, $`\kappa `$ is a propagation constant parameterizing the family of solutions sought for, and $`s`$ is an integer spin. The real amplitude $`U`$ obeys the equation

$$\left(U_{rr}+r^1U_rs^2r^2U+U_{TT}\right)\kappa U+U^3U^5=0,$$
(3)
supplemented by the condition that $`U`$ must decay exponentially as $`r\mathrm{}`$ and $`T\mathrm{}`$ (due to the definition of $`\tau `$ in Eq. (1), the latter condition implies that the solution’s snapshot taken at $`t=\mathrm{const}`$ is localized in the propagation coordinate $`z`$ \- in fact, exactly the same way as 1D solitons are localized in optical fibers Agr ).

Equation (2) conserves the energy $`E=\left|u(X,Y,T)\right|^2𝑑X𝑑Y𝑑T`$, Hamiltonian

$$H=\left[|u_X|^2+|u_Y|^2+|u_T|^2(1/2)|u|^4+(1/3)|u|^6\right]𝑑X𝑑Y𝑑T,$$
(4)
momentum (equal to zero for the solutions considered), and angular momentum in the transverse plane, $`L=\left(\varphi /\theta \right)|u|^2𝑑X𝑑Y𝑑T`$, where $`\varphi `$ is the phase of the complex field $`u`$. Relations between $`L`$, $`H`$ and $`E`$ for a stationary spinning soliton follow from Eq. (3): $`L=sE`$; $`H=\kappa E\frac{2}{3}2\pi rU^6(r,T)𝑑r𝑑T`$ new2 .

We have numerically found families of 3D spinning solitons with a toroidal shape. To quantify the solutions, in Fig. 1 we show the propagation constant $`\kappa `$ and $`H`$ vs. $`E`$ for both $`s=0`$ and $`s=1,2`$ solitons. They exist for $`E`$ exceeding a threshold value, which increases with $`s`$. The full and dashed lines in Fig. 1 correspond to stable and unstable branches according to results presented below. The $`s=0`$ branch of the solutions is divided into stable and unstable portions on the basis of the known criterion which states that the fundamental ($`s=0`$) soliton branch undergoes a change in the stability where $`dE/d\kappa =0`$ VK . However, this criterion ignores azimuthal instability, which is frequently fatal for spinning solitons.

The most revealing information on the stability of solitons is provided by analysis of a linearized version of Eq. (2). To this end, we seek perturbation eigenmodes of the general form,

$`u(Z,r,T,\theta )U(r,T)\mathrm{exp}\left[i(s\theta +\kappa Z)\right]`$ $`=`$ $`f(r,T)\mathrm{exp}\left\{\lambda _nZ+i[(s+n)\theta +\kappa Z]\right\}`$ (5)
$`+g^{}(r,T)\mathrm{exp}\left\{\lambda _n^{}Z+i[(sn)\theta +\kappa Z]\right\},`$
where $`n>0`$ is an arbitrary integer azimuthal index of the perturbation, $`\lambda _n`$ is the (complex) instability growth rate sought for, and the functions $`f`$ and $`g`$ obey equations

$`i\lambda _nf+{\displaystyle \frac{^2f}{T^2}}+{\displaystyle \frac{^2f}{r^2}}+r^1{\displaystyle \frac{f}{r}}(s+n)^2r^2f\kappa f+\left(23U^2\right)U^2f+\left(12U^2\right)U^2g`$ $`=`$ $`0,`$
$`i\lambda _ng+{\displaystyle \frac{^2g}{T^2}}+{\displaystyle \frac{^2g}{r^2}}+r^1{\displaystyle \frac{g}{r}}(sn)^2r^2g\kappa g+\left(23U^2\right)U^2g+\left(12U^2\right)U^2f`$ $`=`$ $`0.`$ (6)
The solutions must decay exponentially at $`r\mathrm{}`$, and vanish as $`r^{\left|s\pm n\right|}`$ at $`r0`$.

To solve Eqs. (6), we used a known numerical procedure unstable1 ; Akh , which produces results presented in Fig. 2. The most persistent unstable eigenmode has $`n=2`$, for both $`s=1`$ and $`s=2`$. As is seen in Fig. 2, with the increase of $`\kappa `$, the instability of the soliton with $`s=1`$, accounted for by $`\mathrm{Re}\lambda _2`$, disappears at $`\kappa =\kappa _{\mathrm{st}}0.13`$, and the stability region extends up to $`\kappa =\kappa _{\mathrm{offset}}^{(3\mathrm{D})}0.17`$, corresponding to infinitely broad solitons (which implies that the vortex of the dark-soliton type vortex , that may be regarded as an infinitely broad spinning soliton, is stable too). The relative width of the stability region is $`\left(\kappa _{\mathrm{offset}}^{(3\mathrm{D})}\kappa _{\mathrm{st}}\right)/\kappa _{\mathrm{offset}}^{(3\mathrm{D})}0.2`$. However, there is no stability region for 3D solitons with $`s=2`$, in contrast to the 2D vortex solitons in the CQ model IsaacPLA . In the case when a spinning soliton is unstable, its instability is oscillatory; the corresponding frequency, $`\mathrm{Im}\lambda `$, is of the same order of magnitude as $`\mathrm{Re}\lambda `$ at the maximum-instability point (see Fig. 2), and $`\lambda `$ becomes purely imaginary at $`\kappa =\kappa _{\mathrm{st}}`$.

The above results were checked in direct simulations of Eq. (2) by means of the Crank-Nicholson scheme combined with the Gauss-Seidel iteration procedure. In Fig. 3 we show the amplitude and energy vs. $`Z`$ for the soliton with $`s=1`$, generated by two different initial configurations, with the same energy of the initial configuration, $`E_0=`$ $`13070`$. Robustness of the spinning STS is attested to by the fact that it can be generated from a Gaussian with a nested vortex whose shape is far from the soliton’s exact form. (Energy loss evident in Fig. 3(b) is caused by emission of radiation in the course of the formation of the stable STS; naturally, the loss is larger for the initial Gaussian configuration, which is farther from the exact soliton’s shape.) Figure 4 shows the gray-scale contour plots of the intensity and phase of both the input Gaussian with a nested vortex and emerging spinning STS at $`Z=400`$.

The instability of the $`s=2`$ solitons is illustrated in Figs. 5 and 6. The azimuthal instability breaks them into zero-spin solitons which fly out tangentially, relative to the circular crest of the original soliton. It is noteworthy that, at an early stage of the evolution shown in Fig. 5, the spinning soliton splits into two fragments, in accordance with the growth-rate calculations predicting the dominant instability to be against the $`n=2`$ perturbation mode, but the subsequent nonlinear evolution is more involved. In Fig. 6(b), the spinning STS splits into three fragments with unequal energies. The fragmentation of this soliton (having $`\kappa =0.09`$) into three parts is in accordance with the growth-rate calculation, see Fig. 2(b).

In conclusion, we have found the first example of stable three-dimensional spinning solitons in a dispersive medium which combines cubic and quintic nonlinearities. Only sufficiently broad solitons with spin $`s=1`$ may be stable. However, the existence of stable spinning 3D solitons is a generic fact, as it is not limited to the cubic-quintic nonlinearity: our preliminary studies indicate that these stable physical objects may also occur in media with competing quadratic and self-defocusing cubic nonlinearities. In fact, a condition which turns out to be necessary for the existence of stable spinning solitons is a competition between two different nonlinearities, one focusing and the other one defocusing.

###### Acknowledgements.

D. Mihalache, D. Mazilu and L.-C. Crasovan acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) and European Community (Access to Research Infrastructure Action of the Improving Human Potential Program). B.A.M. and I.T. appreciates support from the Binational (US-Israel) Science Foundation through the grant No. 1999459. I.T. and A.V.B. acknowledge financial support of the Australian Research Council. L.T. and J.P.T. acknowledge support by TIC2000-1010.


text:
# A SAURON study of M32: measuring the intrinsic flattening and the central black hole mass

## 1 Introduction

M32 is a high-surface brightness, compact E3 companion of the Andromeda galaxy. Ground-based kinematic measurements of this galaxy (Tonry 1984, 1987) already showed a steep gradient in the central velocity profile and a central dispersion peak, suggesting the presence of a central compact object, presumably a supermassive black hole. Since then, (spectroscopic) observations with increasing spatial resolution, both ground-based (Dressler & Richstone 1988; van der Marel et al. 1994a; Bender, Kormendy & Dehnen 1996) and with HST (Lauer et al. 1992; van der Marel, de Zeeuw & Rix 1997; Joseph et al. 2001), strengthened the case for this black hole.

Simultaneously, dynamical models of ever improving quality were built to match these data sets (Richstone, Bower & Dressler 1990; van der Marel et al. 1994b; Qian et al. 1995; Dehnen 1995). For most galaxies, the value of the best-fit black hole mass that is found depends on the assumptions that are made about the distribution function (DF) of the galaxy. Two-integral models, which assume that the DF depends only on the two classical integrals of motion (the energy $`E`$ and the vertical component of the angular momentum $`L_z`$), generally tend to overpredict the central black hole mass, since they cannot provide sufficient radial motion (Magorrian et al. 1998; Gebhardt et al. 2000; Bower et al. 2001).

The current state-of-the art models allow for the maximal degree of anisotropy by assuming that the distribution function depends on three integrals of motion (van der Marel et al. 1998, hereafter vdM98; Cretton et al. 1999, hereafter C99; Gebhardt et al. 2000, 2001). The internal kinematical structure of the best-fitting three-integral model of M32 closely resembles that of an $`f(E,L_z)`$ model. This explains why the central black hole mass that was found in M32 by the early models does not differ much from the current value of $`(3.4\pm 0.7)\times 10^6M_{}`$ (vdM98).

This means that the mass of the central black hole and intrinsic kinematical structure of M32 are constrained rather well. The main uncertainty that remains is the inclination along which M32 is observed, or, equivalently, the intrinsic flattening: edge-on models to a data-set consisting of FOS kinematics and ground-based observations along four position angles fit equally well as models that are projected over $`55^{}`$ (vdM98). When we combine this with the observed flattening of M32 (which is almost constant and equal to $`0.73`$ inside the central ten arcseconds, vdM98), we see that the intrinsic flattening of M32 is not very tightly constrained and can have any value between $`0.55`$ and $`0.73`$.

Recent developments in instrument design offer a way out: integral-field spectrographs capture the full two-dimensional behaviour of objects and are therefore expected to better constrain parameters such as the inclination. This tighter constraint on the intrinsic parameters can be explained from the behaviour of two-integral models, even though these generally provide less accurate fits to observational data. The part of a two-integral distribution function that is even in the velocities is only determined once the meridional plane density distribution is known (Qian et al. 1995). The full intrinsic velocity map is needed to constrain the odd part of the DF, which means we need to measure the line-of-sight velocity $`v_{\mathrm{los}}`$. Whereas these two lowest order velocity moments are sufficient to determine the stellar DF of a galaxy, the dark matter distribution is only constrained once the velocity anisotropy is known (Dejonghe 1987; Gerhard 1991, 1993). This quantity can be measured through the higher order velocity moments. It is therefore necessary to determine the full two-dimensional kinematic behaviour of a galaxy. Capturing this information with only a few slit positions, as is generally attempted, is sometimes possible, but not always, which implies that the distribution function remains unconstrained.

These effects are even more important for three-integral distribution functions, which, in most cases, provide better fits to the observations and have more freedom. Therefore, when galaxy models are constrained by kinematic observations along a modest number of slits, the intrinsic structure may remain largely unconstrained. In this paper, we revisit M32 using two-dimensional kinematical maps obtained with the panoramic integral-field spectrograph SAURON (Bacon et al. 2001, hereafter Paper I; de Zeeuw et al. 2002, hereafter Paper II) and high-resolution spectra obtained with STIS on board HST (Joseph et al. 2001). We show that the use of integral-field data places very tight constraints on the central black hole mass and mass-to-light ratio of M32, but also, for the first time, on the inclination, or equivalently, the intrinsic flattening.

The paper is organized as follows: in Section 2, we give a brief summary of the data that are used to constrain the dynamical models, which are described in Section 3. The results are presented in Section 4 and we show in Section 5 that the use of integral-field data is crucial for most of the results that we obtained. The properties of the best-fitting model are described in Section 6, and we conclude with a discussion in Section 7.

Throughout the paper we assume a distance of $`D=0.7`$ Mpc (Welch et al. 1986). This choice does not influence our conclusions, but sets the scale of the models in physical units. Specifically, lengths and masses scale as $`D`$, while mass-to-light ratios scale as $`D^1`$.

## 2 Observations

We use two fully complementary data sets obtained with current state-of-the-art instruments: the integral-field spectrograph SAURON mounted on WHT, and STIS on board HST. Integral-field spectrographs provide the full kinematic picture in one consistent data-set, while the high-resolution instrument STIS is ideal to probe the central kinematics of nearby objects. With this combined data-set, we expect to be able to determine the intrinsic properties of M32 with great accuracy.

### 2.1 SAURON observations

SAURON is mounted on the WHT on La Palma and is specifically designed to measure the kinematics and line-strengths of nearby ellipticals, lenticulars and spiral bulges. Details on the instrument and the data reduction process can be found in Paper I. The instrument has two observing modes: in the low-resolution mode, the field of view is $`30\times 41`$ arcsec at a pixel size of $`0\stackrel{}{.}94`$, while the high-resolution mode, used to zoom in on objects in good seeing conditions, has a pixel size of $`0\stackrel{}{.}28`$, corresponding to a field-of-view of $`9\times 11`$ arcsec.

The SAURON observations were made on October 15, 1999 (Paper II). Relatively good seeing (FWHM$`0\stackrel{}{.}95`$) allowed us to obtain the measurements in the high-resolution mode. The kinematical maps of M32, derived from a single 2700s exposure, are shown in the four upper panels of Fig. 1. The velocity field reaches a peak value of $``$ 65 km s-1; the misalignment between the rotation axis and the projected minor axis is less than two degrees, demonstrating that the measurements are consistent with axisymmetry. The dispersion of M32 drops well below the instrumental resolution ($`\sigma `$ 90 km s-1) outside the inner few arcseconds. This results in a fourth-order Gauss-Hermite moment $`h_4`$ (measuring the symmetric deviations from a Gaussian, van der Marel & Franx 1993) that is difficult to measure. As is apparent from the fourth panel of Fig. 1, it is almost constant over the field.

A comparison of the SAURON measurements with ground-based long-slit data, obtained with the ISIS spectrograph on WHT, shows that the velocities differ only by $`2.3\pm 2.2`$ km s-1, the dispersions by $`6.0\pm 2.6`$ km s-1 and $`h_3`$ by $`0.009\pm 0.016`$ (Paper II; the accuracy of $`h_4`$ is difficult to determine). This indicates that there is no systematic offset in $`V`$ and $`h_3`$, while the disagreement between the dispersion measurements is only significant at the 2$`\sigma `$-level. We do not correct for this and use all SAURON measurements to constrain the models.

### 2.2 STIS observations

The STIS measurements (Fig. 2) were obtained by Joseph et al. (2001), using the FCQ-method (Bender 1990). They clearly show the main signatures of a central black hole: a steep velocity gradient and a peak in the dispersion profile near the center. By assuming that M32 has a two-integral distribution function and by solving the Jeans equations, Joseph et al. estimate that a central black hole with a mass of $`(24)\times 10^6M_{}`$ best reproduces these signatures, in full harmony with the earlier results by vdM98. STIS provides a continuous spatial sampling at $`0\stackrel{}{.}051`$ pixels, avoiding positioning problems that may occur when using aperture spectrographs such as FOS (van der Marel, de Zeeuw & Rix 1997). Comparisons between these two types of data are therefore only possible if such uncertainties are taken into account (Joseph et al. 2001). The STIS measurements are superior to the FOS data, since the spectral resolution and signal-to-noise ratio are much higher, an advantage when measuring the relatively low velocity and dispersion of M32. Because of this, we can use up to the fourth order Gauss-Hermite moment, while FOS was only able to measure the first two.

## 3 Method

### 3.1 Mass models

We construct dynamical models, based on Schwarzschild’s orbit superposition method (Schwarzschild 1979; 1982). We assume that the galaxy is axisymmetric and allow for the full degree of anisotropy by assuming that the distribution function depends on three integrals of motion. The implementation that is used here was originally developed by Rix et al. (1997), vdM98 and C99, and was designed for use on one-dimensional surface-brightness profiles (but see Cretton & van den Bosch 1999). However, for M32, as for most nearby ellipticals, high-quality two-dimensional imaging is available, so that more sophisticated density parametrisations are possible. An example of such a formalism is Multi-Gaussian Expansion (MGE, Monnet, Bacon & Emsellem 1992; Emsellem, Monnet & Bacon 1994), in which the set of Gaussians that best fits an image of the galaxy is found. The advantages of this formalism are that both axisymmetric and triaxial galaxies can be reproduced in a straightforward manner, the deprojection is simple and can be done analytically (under certain assumptions) once the viewing angles are known, and convolution with (MGE) point spread functions is fast and simple. We have therefore adapted the axisymmetric Schwarzschild software to use MGE; the new algorithm was tested while modeling the dynamics of the counter-rotating core galaxy IC1459 (Cappellari et al. 2002). These tests have shown that the modeling results do not depend on the choice of mass model: the best-fitting parameters and intrinsic (phase space) structure of models that start from an MGE parameterisation with constant axial ratio are identical to those of models that use a power-law density profile.

An MGE parameterisation of M32 was obtained using the software of Cappellari (2002), by fitting to a combination of a high-resolution HST/WFPC2/F814W image and a ground-based $`I`$-band image obtained by Peletier (1993; $`0\stackrel{}{.}549`$ pixels-1 with a seeing of $`1\stackrel{}{.}27`$) with the Isaac Newton Telescope on La Palma. Since we are constructing axisymmetric models, we do not allow for shifts in the position angles of the Gaussians and fix the centers of all individual Gaussians to the center of the galaxy. Each Gaussian therefore has three free parameters: the amplitude $`P^{}`$ (in units of $`L_{}\mathrm{pc}^2`$), the width $`\sigma ^{}`$ in pixels (or arcseconds) and the projected flattening $`q^{}`$. Eleven Gaussians were used in the fit, since, upon adding more Gaussians, the RMS error of the fit does not change by more than 1%. The parameters of these eleven Gaussians are listed in Table 1; this superposition reproduces the M32 surface brightness with an RMS error of $``$2.2 per cent (Cappellari 2002, Figure 6). This is within the measurement errors.

### 3.2 Axisymmetric three-integral Schwarzschild models

We use the MGE-parameterisation to calculate axisymmetric three-integral models with the modified Schwarzschild method. Our models are different from those of previous authors: we use independent kinematic data sets and we parameterise the mass density differently. These differences provide a useful test: if we find best-fitting parameters that match the ones found by e.g. vdM98, it means that both approaches are reliable and the results are robust.

The set of best-fitting Gaussians is deprojected by assuming that the galaxy is axisymmetric and choosing a value of the inclination111 The flattest Gaussian in the MGE superposition, in our case the eighth Gaussian, limits the range of possible inclinations to $`i47^{}`$. Since this range includes values of $`i`$ that can be ruled out for other physical reasons (see below), this mathematical limitation is not a problem.. The resulting intrinsic luminosity density is multiplied by a constant mass-to-light ratio $`M/L`$ to obtain the intrinsic mass density. The stellar potential can then be calculated by applying Poisson’s equation. A central supermassive black hole is mimicked by adding the potential of a central point mass $`M_{}`$. In the combined stellar and black hole potential, a representative orbit library is found by sampling the two classical integrals of motion, the energy $`E`$ and the vertical component of the angular momentum $`L_z`$, together with the third integral $`I_3`$.

Here, we only give a brief outline of the scheme that is used to sample these integrals; a more thorough description can be found in C99. The orbital energy is sampled through a logarithmic grid in the radius $`R_c`$ of the circular orbit at energy $`E`$. This grid is chosen to include more than $`99.99`$ per cent of the mass. While some of the orbits that make up the missing mass might contribute to the density in our radial range, this contribution will be smaller than $`0.01`$ per cent and can therefore be safely ignored.

At each energy, $`L_z`$ is sampled by a linear grid in the fractional angular momentum $`L_z/L_{\mathrm{max}}1,1`$, with $`L_{\mathrm{max}}`$ the maximum orbital angular momentum at energy $`E`$. This grid does not include $`L_z/L_{\mathrm{max}}=\{0,1\}`$ to avoid problems with the numerical integration. This is not a severe limitation, as the orbits that are missed in this manner carry only a small fraction of the mass in an axisymmetric model. The third integral is parameterised by the orbital starting point on the zero velocity curve.

Every combination of the integrals $`(E,L_z,I_3)`$ defines a separate stellar orbit. Since the third integral is generally not known analytically, the orbital properties can only be obtained by solving the equations of motion numerically. At every time step of this orbit integration, which is carried out for a fixed number of characteristic orbital periods (usually chosen to be the period of the circular orbit at the given energy), the orbital energy is randomized by drawing a number of arbitrary values from the energy bin around $`E`$, which are then assigned to the orbit. This ’dithering’ of integral space results in smoother orbital superpositions, which generally provide better fits to the data (Zhao 1996; C99).

The orbital superposition can be made even smoother by including the so-called two-integral components (TICs, C99; Verolme & de Zeeuw 2002) in the orbit library. These TICs are building blocks that obey only the two classical integrals of motion $`E`$ and $`L_z`$ and implicitly include stochastic orbits at the given values of energy and angular momentum. Since they are smoother than most orbits, the orbital superposition also suffers less from discretisation when they are included in the library. Here, we calculated model fits to different data-sets (both long-slit and integral-field data) with and without TICs in the orbit library, and found that the fit improves only marginally upon inclusion of the TICs (RMS error of the fit changes by less than 0.1%). On the other hand, calculating observables for the TICs at every aperture can be very time-consuming when many data points are included in the fit, as is the case when using integral-field data. To save computation time, we therefore decided not to include the TICs in the derivation of the best-fit parameters.

The observables of all orbits (and TICs) in the library are stored on grids that are adapted to the resolution of the observations. Using these, the orbital superposition that best fits the data is found by applying the non-negative least-squares (NNLS) routine of Lawson & Hanson (1974). The (unavoidable) discrepancy between model predictions and data, caused by noise in the data, discretisation effects and wrong choices of the model parameters, is expressed in a value for the $`\chi ^2`$, defined as

$`\chi ^2={\displaystyle \underset{i=1}{\overset{N_d}{}}}\left({\displaystyle \frac{D_i^{}D_i}{\mathrm{\Delta }D_i}}\right)^2,`$ (1)
where $`N_d`$ is the number of constraint points, $`D_i`$ is the observational constraint at the $`i`$-th data point, $`D_i^{}`$ is the model prediction at that point and $`\mathrm{\Delta }D_i`$ is the uncertainty that is associated with this value (usually the observational error). By varying the inclination, central black hole mass and stellar mass-to-light ratio, we investigate which combination $`(M_{},M/L,i)`$ results in the overall smallest $`\chi ^2`$ (the overall best fit to the data). The value of $`\chi ^2`$ for a single model is of limited value, since the true number of degrees of freedom is generally not known, but the difference in $`\chi ^2`$ between a model and the overall minimum value, $`\mathrm{\Delta }\chi ^2=\chi ^2\chi _{\mathrm{min}}^2`$, is statistically meaningful (see Press et al. 1992). By using the central limit theorem, it can be shown that this interpretation is even valid when the errors $`\mathrm{\Delta }D_i`$ are not Gaussian distributed (van der Marel et al. 2000), which is important when applying the statistic to models of observational data. This means that we can assign the usual confidence levels to the $`\mathrm{\Delta }\chi ^2`$ distribution and determine the probability that a given set of model parameters will occur.

Introducing integral-field data in numerical models also means a considerable increase in the number of measurements that has to be reproduced: even when only one SAURON pointing is used, a fit to both the intrinsic and projected mass density and to all four Gauss-Hermite moments results in approximately 8000 constraints. Standard mathematical practice in minimization problems is to use more free parameters than constraints, which means we would need $`10^4`$ orbits to fit a SAURON field. On the other hand, the behaviour of models that are only constrained by observational data is very similar to that of models with additional regularisation constraints (R97; C99; vdM98; Fig. 5 of this paper), even though regularised models always have more constraints than orbits (§4.2).

We conclude from this that it is not necessary to use more orbits than constraints to obtain accurate and unbiased fits to the data. In our case, we decided to use an orbit library with $`N_o=20\times 7\times 7`$ orbits. Since the orbital $`L_z`$ can be both positive or negative, depending on the chosen direction of rotation around the symmetry axis, this results in a library of 1960 orbits, implying that we have almost four constraints per orbit.

When using integral-field data, the integration time of individual orbits may become extremely lengthy, as their properties have to be stored in many apertures. This is another reason not to use very large orbit libraries. Integrating all orbits in our library requires about eight hours of computation time on a 1GHz machine with 512Mb memory, while one NNLS fit takes approximately one hour on the same machine.

## 4 Results

### 4.1 Best-fit intrinsic parameters

We sample the model parameters $`(M_{},M/L,i)`$ on a grid of 990 models, defined by fifteen values of the central black hole mass, spaced linearly between $`M_{}=0`$ and $`M_{}=4.5\times 10^6M_{}`$, eleven values of $`M/L`$ between 1.5 and 2.5, and inclinations $`i=90^{},75^{},70^{},65^{},60^{}`$ and $`50^{}`$. We do not consider smaller values of the inclination, as this would make M32 intrinsically flatter than E7, i.e., flatter than any observed elliptical galaxy. Models with equal ratios between the mass-to-light ratio and the central black hole mass differ only by a scaling factor in the velocity. This means that, for every inclination, only fifteen orbit libraries have to be calculated, each with a different ratio $`M_{}/(M/L)`$. Any point along a line of constant $`M_{}/(M/L)`$ can then be obtained by scaling the velocity prior to the NNLS fit. As a result, one full $`\chi ^2(M_{},M/L)`$ grid (at fixed inclination) takes about twelve days of computation time.

For each of the 990 models that were obtained in this manner, we calculated the fit to the combined SAURON and STIS data. Contour plots of the values of $`\chi ^2(M_{},M/L,i)`$, offset such that the overall minimum value is zero, are shown in Fig. 3.

The contour levels correspond to the confidence levels of a $`\mathrm{\Delta }\chi ^2`$ distribution with three degrees of freedom (with the thick contour in the figure indicating the 3$`\sigma `$ confidence level), so that the volume in which the intrinsic parameters of M32 are most likely located can be readily assessed. The shape of the contours is very regular, even though no smoothing has been applied.

The best-fitting edge-on model has a black hole mass of $`(4.0\pm 0.5)\times 10^6M_{}`$ (3$`\sigma `$-error), which agrees very well with the vdM98’s value of $`(3.4\pm 0.7)\times 10^6M_{}`$, even though the observational data sets are independent. However, we see from Fig. 3 that models with an inclination in the range $`i=70^{}\pm 5^{}`$ provide significantly better fits to the data than models with other inclinations, which are ruled out at high confidence. This narrow constraint on the inclination corresponds to an intrinsic flattening of approximately $`0.68\pm 0.03`$ (given the observed flattening of 0.73).

The overall best-fit model, located at $`i=70^{}`$, has a black hole mass of $`(2.5\pm 0.5)\times 10^6M_{}`$ and an $`I`$-band mass-to-light ratio of $`(1.85\pm 0.15)M_{}/L_{}`$. The fit of this model to the SAURON data is shown in the bottom panels of Fig. 1; the fit to the STIS data is given in Fig. 2.

### 4.2 Applying regularisation

The matrix problem that is solved by the NNLS routine is often ill-conditioned, resulting in orbital weights $`\gamma (E,L_z,I_3)`$ that vary rapidly as a function of the integrals of motion. This can be avoided by imposing regularisation constraints on the orbital weights (see e.g. Press et al. 1992; Merritt 1993; Zhao 1996; Rix et al. 1997); another method, with similar effect, is to add maximum entropy constraints (Richstone & Tremaine 1988). Applying regularisation forces the orbital weights towards a smooth function of $`E,L_z`$ and $`I_3`$ by minimizing the $`n`$-th order derivatives $`^n\gamma (E,L_z,I_3)/E^n,^n\gamma (E,L_z,I_3)/L_z^n`$ and $`^n\gamma (E,L_z,I_3)/I_3^n`$. The degree of smoothing is determined by the order $`n`$ and by the maximum fractional error $`\mathrm{\Delta }`$ that the derivatives are allowed to have. For all regularised models that are described here, we have minimized the second-order derivative ($`n=2`$).

The optimal value of the regularisation error $`\mathrm{\Delta }`$ was determined by calculating model fits for a fixed combination of ($`M_{},M/L,i`$) and varying $`\mathrm{\Delta }`$ between $`0.1`$ (high regularisation) and $`\mathrm{}`$ (no regularisation). Fig. 4 shows that $`\chi ^2`$ is large for small values of $`\mathrm{\Delta }`$, but quickly decreases to an almost constant level for $`\mathrm{\Delta }>4`$.

Since the ideal model both fits the constraints and has a smooth distribution of orbital weights, we use $`\mathrm{\Delta }=4`$ whenever regularisation is needed. A similar degree of smoothing was used earlier in applications that use the same regularisation technique (vdM98; Cretton & van den Bosch 1999) and in models that use the maximum entropy method (Gebhardt et al. 2000).

Including smoothing constraints in the NNLS fit will, in general, change the value of $`\chi ^2`$ and therefore may affect the shape of the contours in plots such as Fig. 3. The top panel of Fig. 5 shows a cross-section of Fig. 3 at $`M/L=1.8`$; the $`\chi ^2`$ of fits with the same parameters, but including regularisation, is shown in the middle panel. Since the shape of the contours does not change significantly, applying regularisation does not change the overall best-fit parameters, while it allows us to better study the phase-space structure and internal motions (see Fig. 7 below).

The two upper panels were calculated using four constraints per orbit ($`N_d/N_o=4`$, the value that is used throughout the paper). Although we already argued in § 3.2 that our results are not sensitive to this parameter, we investigated this more quantitatively by calculating (unregularised) models for which $`N_d/N_o=0.6`$, which means there are almost two orbits per constraint point. The result of this is shown in the bottom panel, which indeed illustrates that our results do not depend on the number of orbits that is used.

## 5 The effect of integral-field data

At this point, although we have shown that the model parameters are constrained much more tightly than in any previous model, it is not entirely clear to what extent this is caused by the use of integral-field data. We also have to take into account that we have used a very dense sampling of $`M_{}`$ and $`M/L`$, that we were able to calculate models at more inclinations than vdM98, and that the STIS data is of higher quality than the FOS observations. A more quantitative estimate of the influence of the integral-field data on the $`\chi ^2`$ contours and on the best fitting parameters can be obtained by repeating the procedure of Section 3 for a data set that closely resembles the one used by vdM98. Since they fitted to a combination of FOS pointings with ground-based long-slit observations along the major, minor, $`45^{}`$ and $`135^{}`$ axes (which agree with SAURON), a similar data set can be obtained by extracting virtual ‘slits’ from the SAURON field along the same position angles and combining these with the STIS data.

This data set was used to constrain models with the same parameter ranges as described in Section 3. Fig. 6 shows the $`\chi ^2`$ of these fits. The best-fitting black hole mass is $`M_{}=(2.8\pm 1.2)\times 10^6M_{}`$, the $`I`$-band mass-to-light ratio is $`(1.75\pm 0.25)M_{}/L_{}`$ and the inclination can have any value below $`75^{}`$. A comparison with the best-fit parameters of Section 4 indicates that the volume of intrinsic parameters providing a reasonable fit to the data is significantly larger for the slit data set than when the entire SAURON field is used. We also notice that, although the use of integral-field data tightens the constraint on all model parameters (with a fractional change in the allowed range of $`M_{}`$ of 2.4 when changing from integral-field to slit data, and 1.7 for $`M/L`$), the most dramatic effect occurs for the inclination, which has an eight times larger allowed range when the galaxy is observed along a few slits only.

This can be explained from the behaviour of two-integral models, even though these generally provide less accurate fits to observational data. As was mentioned in the Introduction, in order to constrain the even part of a two-integral distribution function, one needs to specify the meridional plane density distribution (Qian et al. 1995). The odd part is determined once the full intrinsic velocity field is known, which can only be measured through the projected velocity field $`v_z^{}(x^{},y^{})`$ (in which $`x^{}`$ and $`y^{}`$ are the coordinates in the plane of the sky). Furthermore, even if the total stellar DF is determined accurately, the dark matter distribution is still unconstrained. This can be solved by measuring the higher order velocity moments, which are related to the anisotropy of the velocity distribution.

When the projected velocity moments are simple functions of $`x^{}`$ and $`y^{}`$, measurements at a small number of positions are sufficient to determine the behaviour at all intermediate position angles. Indeed, some very specific two-integral models, such as the power-law galaxies (Evans & de Zeeuw 1994), have velocity fields that are quadratic in $`x^{}`$ and $`y^{}`$, implying that observations along two perpendicular slits are sufficient to determine their full two-dimensional kinematical behaviour.

However, most galaxies require a distribution function that depends on three integrals of motion. The velocity fields of these three-integral models (and, in fact, of most generally applicable two-integral models) are usually much more complicated functions of the projected coordinates (Dehnen & Gerhard 1993; de Bruijne, van der Marel & de Zeeuw 1996). This means that these velocity fields cannot be characterized by a (small) number of slit orientations. It is therefore not surprising that realistic models constrained by observations along only a few slit orientations cannot provide much information about the intrinsic shape of the galaxy. It also implies that integral-field data places constraints on the intrinsic shapes of galaxies, as well as on their mass distribution.

We see in both Fig. 3 and Fig. 6 that the best-fitting black hole mass of each separate panel increases with inclination. This behaviour can be understood at least partially from the following: the contribution of the intrinsic rotational velocity $`v_\varphi `$ to the line-of-sight velocity decreases proportionally to $`\mathrm{sin}i`$. However, at smaller inclinations, the model must also be intrinsically flatter, which means $`v_\varphi `$ is larger when lowering $`i`$. For a two-integral model, the increase in $`v_\varphi `$ is dominant, so that there is more net velocity at smaller values of the inclination. This also implies that less mass is needed to obtain the same velocity field. Unfortunately, since there is a trade-off between the mass-to-light ratio and the central black hole mass and because we are dealing with three-integral models, it is not clear to what extent these considerations can explain the effect we observe.

## 6 Properties of the best-fit model

The upper panels of Fig. 7 show the phase-space distribution of the orbits of the non-regularised best-fit model as a function of $`I_3`$ and $`L_z`$, for constant $`E`$. The energy decreases from left to right, corresponding to an increasing distance from the galaxy center. The bottom panels are similar, for the model with a regularisation error of $`\mathrm{\Delta }=4`$.

We can study the intrinsic dynamical structure of the best-fitting model ($`M_{}=2.5\times 10^6M_{}`$, $`M/L`$=1.85 $`M_{}/L_{,I}`$, $`i=70^{}`$) by adding the appropriate moments of the orbits with non-zero weight. This implies that we can study the degree of anisotropy of the model, and from this determine to what extent the galaxy resembles a two-integral model. Fig. 8 shows the intrinsic velocity moments $`v_r^2^{\frac{1}{2}}`$ (solid line), $`v_\theta ^2^{\frac{1}{2}}`$ (dotted line) and $`v_\varphi ^2^{\frac{1}{2}}`$ (dashed line) of the best-fitting model, as a function of the meridional plane radius and the angle from the symmetry axis. The velocity moments are normalized by the total RMS motion $`v^2=v_r^2+v_\theta ^2+v_\varphi ^2`$. We see that the model resembles an $`f(E,L_z)`$ model near the equatorial plane ($`v_r^2^{\frac{1}{2}}=v_\theta ^2^{\frac{1}{2}}`$) and becomes increasingly more anisotropic towards the symmetry axis. This indicates that a two-integral model is not suitable to fit all constraints simultaneously, in agreement with the findings of vdM98 (their Fig. 10). Furthermore, as in their models, azimuthal motion dominates.

For completeness, we have checked the internal kinematical structure of the pseudo-slit model for the same parameters as were used for Fig. 8. Since this combination of intrinsic parameters is inside the 3$`\sigma `$ contour of Fig. 6, and therefore also provides a good fit to the data, the results resemble those of Fig. 8 closely. The main difference between using slits and integral-field data is that more combinations of intrinsic parameters can be ruled out by integral-field data. A point ($`M_{},M/L,i`$) that is not inside the 3$`\sigma `$ contour of Fig. 3, but cannot be ruled out by observations along a few slits, will therefore in general show a different intrinsic kinematical structure.

The anisotropy of a model can also be expressed in terms of an anisotropy parameter. A possible definition of this parameter is

$`\beta _B=1{\displaystyle \frac{v_\theta ^2^{\frac{1}{2}}+v_\varphi ^2^{\frac{1}{2}}}{2v_r^2^{\frac{1}{2}}}}.`$ (2)
When the model is fully isotropic, $`\beta _B=0`$, while $`\beta _B=1`$ for a model that consists entirely of radial orbits, and $`\beta _B=\mathrm{}`$ describes a model with only circular orbits. Fig. 9 shows this parameter as a function of the meridional plane coordinates $`R`$ and $`z`$, over the range that is constrained by the SAURON data. A typical isodensity contour is overplotted (using $`i=70^{}`$). We see that $`\beta _B`$ is positive close to the galaxy center, so that the velocity distribution is radially anisotropic there. In all other parts of the meridional plane, the model is azimuthally anisotropic. This result does not change upon adding regularisation.

Dynamical models of elliptical galaxies have revealed that a few other objects are also tangentially anisotropic (NGC4697, Dejonghe et al. 1996; NGC1700, Statler et al. 1999), while the majority seems to be radially anisotropic (M87, Merritt & Oh 1997; NGC2324, R97; NGC6703, Gerhard et al. 1998; M32, vdM98; NGC1600, Matthias & Gerhard 1999; NGC3379, NGC3377, NGC4473, NGC5845, Gebhardt et al. 2000; NGC1399, Saglia et al. 2000). However, this radial anisotropy is usually rather small and sometimes even accompanied by a transition to isotropy at small radii. Furthermore, not all models included the full velocity profile in the fit, which implies that the mass-anisotropy degeneracy may not have been broken. Finally, our models are the first to use fully two-dimensional kinematical observations, which makes comparison with earlier models even harder. Therefore, although our results seem to contradict these findings, we must also conclude that our knowledge of anisotropy is still limited.

## 7 Conclusions

We have presented dynamical models of the nearby compact elliptical M32, using data from the integral-field spectrograph SAURON and from STIS on board HST. We have shown that our modeling software is able to deal with two-dimensional kinematical information and that the integral-field data tightens the constraint on all intrinsic model parameters considerably.

The axisymmetric three-integral model that best fits the data has a black hole mass of $`(2.5\pm 0.5)\times 10^6M_{}`$ and a stellar $`I`$-band mass-to-light ratio of $`(1.85\pm 0.15)M_{}/L_{}`$. These values confirm the best-fit parameters that were obtained by previous authors, although our modeling procedure is different: we use a fully independent kinematical data-set and a different parameterisation for the mass density. Despite these differences in the methods, the best-fitting parameters that we find match the ones that were found by e.g. vdM98. This means that both approaches are reliable and that the results are robust.

For the first time, we are able to determine the inclination along which M32 is observed with great accuracy: the best-fit value is $`70^{}\pm 5^{}`$. If M32 is indeed axisymmetric, the averaged observed flattening of 0.73 then corresponds to an intrinsic flattening of $`0.68\pm 0.03`$. We have shown in Section 5 that this tight constraint is mainly caused by the use of integral-field data. This implies that integral-field data will provide us with more insight into the internal structure and kinematics of such objects.

Although M32 is consistent with axisymmetry, it may be intrinsically triaxial, but seen along one of its principal planes. Allowing for an intrinsically triaxial object also would enable us to study the effects of uncertainties in the deprojection on the best-fitting parameters more closely. An extension of our method to triaxiality is in progress (Verolme et al. 2002, in prep.), which will allow us to verify the assumptions made here. Furthermore, the tests that were carried out in this paper and the agreement with the results of other authors indicate that, given the assumptions, the results are robust.

We have obtained observations with SAURON of a representative sample of nearby ellipticals, lenticulars and spiral bulges, for most of which high-resolution STIS data is, or will become, available. We are in the process of applying the axisymmetric version of Schwarzschild’s method on SAURON observations of a few of the sample galaxies that appear consistent with axisymmetry (e.g. NGC821, NGC3377, NGC2974). The results of these dynamical models will provide us with the intrinsic parameters of a considerable number of objects, giving us unique insight into the formation and evolution of early-type galaxies.

Acknowledgements
We thank Karl Gebhardt, Richard McDermid and Glenn van de Ven for useful discussions, and Eric Emsellem, Harald Kuntschner and Reynier Peletier for a critical reading of the manuscript. YC acknowledges support through a European Community Marie Curie Fellowship. MB acknowledges support from NASA through Hubble Fellowship grant HST-HF-01136.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.