text:
# The high energy emission line spectrum of NGC 1068

## 1 Introduction

The archetypal Seyfert 2 galaxy, NGC 1068, has been extensively studied at all wavelengths. It was the discovery of broad lines in the polarized flux of this source that led Antonucci & Miller (1985) to propose the Unification model for Seyfert galaxies, which is at present the most popular scenario for this class of sources. In X–rays, after the pioneering $`Einstein`$ and EXOSAT observations (Monier & Halpern 1987; Elvis & Lawrence 1988), the $`GINGA`$ discovery of a strong iron line (Koyama et al. 1989), being interpreted as due to the reprocessing by circumnuclear matter of the otherwise invisible nuclear radiation, brilliantly confirmed the Antonucci & Miller model. ASCA (Ueno et al. 1994; Iwasawa et al. 1997; Bianchi et al. 2001) demonstrated that the line is complex, implying reflection from both neutral and ionized material. The two reflectors model was confirmed, observing the hard X–ray continuum, by BeppoSAX (Matt et al. 1997a), which was also able to put a lower limit of about 1025 cm-2 to the column density of the absorbing matter. Guainazzi et al. (1999) and Bianchi et al. (2001) demonstrated that the situation is even more complex, the line spectrum requiring at least three reflectors, one neutral (cold), one mildly ionized (warm) and one highly ionized (hot). Comparing the two BeppoSAX observations, taken about one year apart, Guainazzi et al. (2000) revealed a flux variability, which they interpreted as due to a variation of the intensity of the hot ionized reflector, so placing an upper limit to its size of the order of a parsec or so.

High spatial and spectral resolution $`Chandra`$ and XMM–$`Newton`$ observations have confirmed the complexity of the circumnuclear region. The $`Chandra`$ image (Young et al. 2001) revealed, especially in soft X–rays, a very rich morphology, with the brightest spot, however, confined in a 1″.5 (corresponding to 118 pc assuming $`H_0`$=70 km/s/Mpc, Young et al. 2001) region around the nucleus. Indeed, the iron emission as well as the neutral and ionized continuum reflection are mostly confined in the nuclear region (Ogle et al. 2003). Many off–center point–like sources are also detected (Smith & Wilson 2003).

Gratings observations (Kinkhabwala et al. 2002; Ogle et al. 2003) have measured a line spectrum which is consistent with photoionized plasma with a wide range of ionization parameters. The line fluxes require significant contribution from resonant scattering (Kinkhabwala et al. 2002), as predicted by Band et al. (1990) and Matt et al. (1996). Given the limited energy range of the XMM–$`Newton`$/RGS and the modest effective area of the $`Chandra`$/HETG at high energies, these results mainly concern the soft X–ray line spectrum.

In this paper we analyse and discuss the high energy ($`>`$4 keV) XMM–$`Newton`$ spectrum of NGC 1068. While the CCDs of the EPIC cameras have a significantly poorer energy resolution than the $`Chandra`$ HETG (the only grating instrument currently working at the iron K$`\alpha `$ energy), for many purposes the much larger effective area above 4 keV overcompensate for the lower resolution. As an example, one can compare the wealth of informations on Fe and Ni K lines derived from the XMM–$`Newton`$ observation of the Circinus Galaxy (Molendi et al. 2003) with the much poorer informations on the same lines (indeed many of them not even detected) derived from the $`Chandra`$/HETG spectrum (Sambruna et al. 2001).

## 2 Data reduction

### 2.1 XMM–Newton

NGC 1068 was observed twice by XMM-Newton on July 29, 2000 and July 30, 2000 with the imaging CCD cameras, the EPIC-pn and MOS, adopting the Medium filter and operating in Large Window (pn), Full Window (MOS1) and Prime Partial W2 (MOS2). We retrieved the ODF files from the public archive, and reduced them with the latest available SAS version (5.4.1), therefore using the most updated calibrations. The observed count rate is lower than the maximum defined for a 1% pileup for the selected pn subarray (see Table 3 of XMM-Newton Users’ Handbook), but it is close to that value, so we conservatively decided to exclude patterns higher than 0 in data reduction. An analysis with the epatplot tool indicates that same residual pile–up effect may still be present, so we extracted a spectrum removing the inner 6″ circle region. No significant differences are found, however, in the continuum and in the lines alike, if the comparison is limited to the energy range of interest here, i.e. above 4 keV. We therefore concluded that, for the pn, pile–up is not relevant for our purposes. On the other hand, we shall not use MOS data because MOS1 has been used in Full Window mode and therefore strongly suffers from pile–up problems, while MOS2 does not add much information, expecially because we will limit our analysis to $`E>`$4 keV. Moreover, the pn appears to be the best calibrated detector at the energies under consideration, as deduced from the analysis of the iron line complex in the spectrum of the Circinus Galaxy (Molendi et al. 2003).

Since the two observations are nearly contiguous and no significant flux variations are found in their lightcurve, we verified that the two spectra were consistent with each other and decided to combine them. We used the SAS tool merge to obtain a single event file and then extracted a spectrum from a region of $`40\mathrm{}`$ radii. The total exposure time, after rejecting time intervals of flaring particle background, is 61 ks. Spectra were analysed with Xspec 11.2.0.

### 2.2 Chandra

Given the complexity of the morphology of the circumnuclear region of NGC 1068, we also analysed the Chandra observation of NGC 1068 performed on February 21, 2000 with the Advanced CCD Imaging Spectrometer (ACIS-S), for a total exposure time of 48 ks. The default 3.2 s frame time, though resulting in the nuclear spectrum spoiled by pileup, does not affect our qualitative imaging analysis and the spectrum of the much fainter source CXJ024238.9-000055 (see below). Data were reduced with CIAO 2.3, using CALDB 2.21. Spectra were analysed with Xspec 11.2.0.

## 3 Spatial analysis

In Fig. 1 the $`Chandra`$ image over the whole working band (left panel) and above 4 keV only (right panel) is shown. The detailed spatial analysis is discussed by Young et al. (2001), while the X-ray source population by Smith & Wilson (2003). The XMM–$`Newton`$ image (whole band, left panel, and above 4 keV, right panel) is shown in Fig. 2. In all images, the $`r`$=40″ circle over wich the XMM–$`Newton`$ spectrum has been extracted is shown.

For our own purposes, it is important to note that: a) the nucleus is partly extended over several arcseconds even above 4 keV, the extended emission accounting for about a quarter of counts in this band. Our 40″ XMM–$`Newton`$ extraction regions encompasses this extended emission. b) At a distance of about 28″, a bright point–like source (CXJ024238.9-000055) is clearly apparent in the $`Chandra`$ image (hereinafter we refer to it for simplicity as X–1); in the hard band, it is by far the brightest off–nuclear point–like source. Smith & Wilson (2003) derived for this source a 0.4-8 keV flux of about 5$`\times 10^{13}`$ erg cm-2 s-1, corresponding to a luminosity, at the distance of NGC 1068, of 1.4$`\times 10^{40}`$ erg s-1. The spectrum (which we reanalysed) is quite hard, being consistent with either a very flat ($`\mathrm{\Gamma }`$=0.85) power law or with multicolor disc blackbody with an inner temperature of about 4 keV. While we defer the reader to Smith & Wilson (2003) for a discussion of its possible nature, we note here that this source, being so hard, may significantly contaminate the spectrum of the NGC 1068 nucleus in the energy range under consideration in this paper. Indeed, while the source is not readily visible in the whole band XMM–$`Newton`$ image, it is clearly present in the $`E>`$4 keV image. To evaluate its contribution to the overall spectrum, we extracted a 20″ spectrum centered on NGC 1068 and a 10″ spectrum centered on X–1. The flux we derived for NGC 1068 is 14 times larger than that of X–1 (which in turn is about half that measured by $`Chandra`$, Smith & Wilson 2003), and this is certainly a lower limit as the spectrum of X–1 is heavily contaminated by that of NGC 1068. In fact, in the X–1 spectrum a broad iron line is clearly present, well fitted by two narrow components, one at 6.4 keV (neutral iron) and one at 6.7 keV (He–like iron). If we assume that they are the result of the contamination by NGC 1068 (no iron lines are apparent in the, admittedly poor, $`Chandra`$ spectrum), this can provide a method to estimate the level of contamination of X-1 by the nucleus. The EWs for these two lines are about 2.5 smaller than in the NGC 1068 spectrum (when calculated with respect to the whole continuum; see next section). Therefore, about 40% of the spectrum extracted around X–1 is actually due to NGC 1068, which implies that X-1 does not contaminate the 40″ spectrum of NGC 1068 by more than about 5%. Considering also that we are mostly interested in the emission lines, which are unlikely to be produced by X–1, we decided that this contamination is not large enough to justify the 20% loss in photons when using the $`r=20\mathrm{}`$ region. In the following we will then use the 40″ radius spectrum, unless explicitely stated.

## 4 The E$`>`$4 keV XMM–$`Newton`$ spectrum

The soft X–ray high energy resolution line spectrum of NGC 1068 is discussed in detail by Kinkhabwala et al. (2002) and Ogle et al. (2003). We limit ourselves to analyse and discuss the high energy ($`E>`$4 keV), moderate resolution EPIC–pn spectrum of the nuclear region of NGC 1068, with particular emphasis on the iron and nickel line properties. Let us, however, first briefly discuss the continuum properties.

### 4.1 The continuum

Following Iwasawa et al. (1997) and Matt et al. (1997a), we fitted the continuum with two components, describing reflection from both neutral and ionized matter (the latter component including both the ‘warm’ and ‘hot’ regions, Bianchi et al. 2001). For the neutral reflection component we adopted, as customary, the model pexrav, with $`R`$ fixed to –1 (i.e. with the primary emission shut off). The ionized reflector component is instead modeled by a simple power law, as the ionized reflector should act as a “mirror” for the primary emission. The power law index has been assumed to be the same for the two models. To fit the spectrum, we of course added as many lines as required, as discussed in details in the next paragraphs. The fit is fully acceptable ($`\chi ^2`$/d.o.f.=143.7/149). The best fit parameters for the continua, for the neutral and the ionized lines are reported in Tables 1, 2 and 3, respectively (all errors refer to 90% confidence level for one interesting parameter, unless explicitely stated). The folded spectrum and best fit model, along with the model/data ratio, are shown in Fig. 3, while the unfolded spectrum is shown in Fig. 4. The power law index is 2.04$`{}_{0.03}{}^{}{}_{}{}^{+0.04}`$, and the ratio between the 4–10 keV cold and ionized reflectors fluxes is about 1.5. (Fitting the spectrum extracted from a 20″ radius, a slightly steeper index, 2.10$`{}_{0.04}{}^{}{}_{}{}^{+0.04}`$, is found; the ratio between the two components does not change significantly. All other parameters are also consistent with those obtained with the $`r`$=40″ extraction radius, but with of course somewhat larger error bars). This result could be compared with that derived by Matt et al. (1997a) from the December 1996 BeppoSAX observation, where the flux ratio was instead 0.5, suggesting significant variations in likely both components. Indeed, variability on a time scale of 1 year, probably in the ionized reflector component, has been reported by Guainazzi et al. (2000). However, the BeppoSAX result must be taken with caution: the $`F_{\mathrm{cold}}`$/$`F_{\mathrm{ion}}`$ ratio is very sensitive to the power law index, which in BeppoSAX was only loosely determined (1.74$`{}_{0.56}{}^{}{}_{}{}^{+0.25}`$). We re–fitted the BeppoSAX spectrum imposing the power law index found from XMM–$`Newton`$, i.e. 2.04. The result is only slightly worse statistically $`\mathrm{\Delta }\chi ^2`$=2.7 with one less free parameter, corresponding to a 90% confidence level according to the F-test), and the ratio now becomes $``$1.3. On the contrary, fitting the XMM–$`Newton`$ spectrum imposing the BeppoSAX best fit power law index gives a significantly worst fit ($`\chi ^2`$/d.o.f.=174.4/150); the ratio is now about 1.15. Indeed, the XMM–$`Newton`$ 90% confidence level lower limit on $`\mathrm{\Gamma }`$ (i.e. 2.01) corresponds to a flux ratio of 1.49. At this confidence level, therefore, the BeppoSAX and XMM–$`Newton`$ estimates are inconsistent with each other. There is therefore some, even if still not conclusive given also the uncertainties when comparing different instruments, evidence for variability between the BeppoSAX and the XMM–$`Newton`$ observations, taken 3.5 years apart, which would give an upper limit to the distance of the spatially unresolved reflecting regions of the order of a few parsecs. Clearly, a second XMM–$`Newton`$ observation is necessary to confirm this finding.

Another interesting result is the large iron abundance, $`A_{\mathrm{Fe}}`$=2.4 (with respect to the solar value of Anders & Grevesse 1989). In the fit, $`A_{\mathrm{Fe}}`$ is basically derived from the depth of the iron edge in the Compton reflection component, which depends on the spectral index and, for a plane–parallel geometry, also on the inclination angle. Solar abundances for the other elements (in particular for CNO elements, which contribute most to the photoelectric cross section at 6.4 keV) have been assumed. The fit has been performed keeping the inclination angle $`i`$ fixed to 63, i.e. the Xspec default value. For lower values of the inclination angle a less deep edge is expected, and therefore an even larger value of the iron abundance is found (e.g. 3.3 for an inclination of about 20). In Figs. 5 the $`A_{\mathrm{Fe}}`$$`\mathrm{\Gamma }`$ (with $`i`$ free to vary) is shown; it is important to note that the iron abundance is larger than 1.5 at the 99% confidence level (two interesting parameters).

### 4.2 Emission lines from the cold reflector

Of the several lines observed in the high energy spectrum of NGC 1068 (see Fig. 3), four can be attributed to neutral Fe and Ni atoms, i.e. the K$`\alpha `$ and K$`\beta `$ lines from both elements (see Table 2). All lines, both neutral and ionized, have been fitted with a gaussian profile with $`\sigma `$ fixed to 1 eV (i.e. much less than the energy resolution; no improvement is found leaving $`\sigma `$ free to vary for the brigthest lines). All lines we attribute to the neutral reflector are significant at more than 99.99% confidence level, according to the F-test, apart from the Ni K$`\beta `$ at 8.25 keV which is significant at the 99.90% confidence level.

The iron K$`\alpha `$ line energy (6.424$`\pm `$0.001 keV; note that the statistical error is much lower than the calibration uncertainty, which is of the order of 10 eV111see http://xmm.vilspa.esa.es/docs/documents/CAL-TN-0018-2-1.pdf.) is larger than the value for neutral iron (actually a doublet with energies of 6.391 and 6.405 keV, House 1969, with a 1:2 branching ratio), and would instead suggest Fe xvi; however, the measured K$`\beta `$/K$`\alpha `$ flux ratio for iron (including only the line core for the K$`\alpha `$ line, see below) is 0.20$`{}_{0.03}{}^{}{}_{}{}^{+0.07}`$, slightly larger than the expected value even for neutral value, which is about 0.16 (see Molendi et al. 2003). This value stays almost constant up to Fe ix, then starts decreasing (Kaastra & Mewe 1993). It is therefore unlikely that iron is much more ionized than that (the line energy is also basically constant over this ionization range, House 1969). The discrepancy with the observed line energy may result from a less than perfect energy calibration of the pn, or a real blueshift of the line, which is however hard to explain in the classical torus model for the neutral reflection. It is worth noting that the MOS2 gives exactly the same best fit line energy, and that $`Chandra`$/HETG also measured a blueshift, albeit smaller (6.411$`\pm `$0.007 keV, Ogle et al. 2003). The K$`\beta `$/K$`\alpha `$ ratio for nickel is 0.57$`\pm `$0.35, definitely larger than the expected value, which is similar to that for iron (Kaastra & Mewe 1993). While the measured energy of the line, 8.25$`{}_{0.03}{}^{}{}_{}{}^{+0.10}`$ keV, is perfectly consistent with the expected value of 8.265 keV for the neutral Ni K$`\beta `$ line (Bearden 1967), the Fe xxvi K$`\beta `$ lines at energies clustering around 8.25 222Unless explicitely stated, line energies are taken from the NIST Atomic Spectra Database, http://physics.nist.gov/cgi-bin/AtData/main\_asd. are also possibly contributing.

The Ni–to–Fe K$`\alpha `$ flux ratio is 0.13$`\pm `$0.03, much larger than the expected value of 0.04-0.06 (calculated as described in Molendi et al. 2003 assuming a factor 2 overabundance for both iron and nickel). This implies a nickel–to–iron overabundance of a factor $``$2, or more

Clear evidence for an iron K$`\alpha `$ Compton shoulder (hereinafter CS) is also found, confirming and refining the ASCA (Iwasawa et al. 1997) and $`Chandra`$/HETG (Ogle et al. 2003) results. The CS is fitted with a gaussian profile with the centroid energy and $`\sigma `$ fixed to 6.3 keV and 40 eV, respectively (Matt 2002). If the CS is not included a significantly worse fit ($`\chi ^2`$/d.o.f.=181.2/149) is found, even allowing for the width of the main line to vary. The ratio between the CS and the line core fluxes is about 0.2, in agreement with the theoretical expectation for Compton–thick reflecting material (Matt 2002), so suggesting an origin in the visible part of the $`N_H>`$1025 cm-2 absorbing material (Matt et al. 1997a).

The iron K$`\alpha `$ EW (with respect to the cold reflection continuum) is 1.2 keV (1.45 keV including the CS), reasonably in agreement with theoretical expectations (e.g. Ghisellini et al. 1994, Matt et al. 2003) even with the observed iron overabundance, given the less than linear increase of the line EW with the iron abundance (Matt et al. 1997b).

### 4.3 Emission lines from the ionized reflector

Four emission lines from highly ionized matter (the ‘hot’ reflector cited above) are also detected in the pn spectrum (Table 3). All of them are significant at more than 99.99% confidence level, according to the F-test, apart from the Ni He–like K$`\alpha `$ at 7.83 keV which is significant at the 99.79% confidence level. Besides the He– and H–like Fe K$`\alpha `$ lines333Also for these lines the measured energies are not in agreement with the expected values, if only statistical errors are taken into account. However, in these cases identification is certain, as no other strong lines are present at those energies. It is worth noting that the energy of the He–like line is greater than expected, and that of the H–like line lower, ruling out inflow or outflow if the two lines, as likely, are produced in the same region., already discovered by ASCA (Ueno et al. 1994; Iwasawa et al. 1997; Bianchi et al. 2001), two lines at about 6.61 and 7.83 keV are also detected with high significance. The energy of the first line points towards fluorescence from Be–like iron, which however has a quite low fluorescent yield of 0.11 (Krolik & Kallman 1987). Moreover, the matter may likely be optically thick to the line, which results in an effective Auger destruction. It is therefore possible that this line is actually a blend of Be– and Li–like lines (the latter having an energy of about 6.65, a high fluorescent yield, 0.75, and no possibility of Auger destruction). (The possibility that this line is the CS of the He–like iron line is ruled out by the low optical depth of the hot reflecting matter, as discussed below). The line at 7.83 keV is likely a blend of the Ni He–like K$`\alpha `$ (7.81 keV) and Fe He–like K$`\beta `$ lines (7.87 keV). The rather wide range of ionization argues against collisional ionization. On the contrary, we have verified with cloudy444http://www.nublado.org/ that photoionization may produce such a range in a single zone with e.g. a column density a few times 1021 cm-2 (see below) and a ionization parameter of about 150.

The measured values of the EWs can be compared with the results of the best previous observation, i.e. the 100 ks ASCA observation (Bianchi et al. 2001). The results reported there are somewhat different, but were obtained with a simpler model. We refitted the ASCA spectrum with all the lines found by XMM–$`Newton`$ (and the power law index fixed to 2.04) and found that, within the errors, the fluxes measured by ASCA are consistent with those measured by XMM–$`Newton`$. It should also be noted that the presence of the 6.61 keV line vindicates the claim of Iwasawa et al. (1997), based on the first ASCA observation, of a redshift of the ionized lines.

The EW of the He– and H–like iron lines can be compared with the results of Matt et al. (1996) to derive at least an order–of–magnitude value for the column density of the reflecting material. In their figures 5 and 6 the EW of these lines, calculated including contributions from both recombination and resonant scattering, are shown as a function of the column density of the reflecting material. The calculations assumed iron abundance and fractions equal to 1 of the considered ions. The measured iron overabundance roughly compensates for the less than one fraction of each ion (for instance, our cloudy calculation gives a value of 0.4 for the fraction of He–like iron atoms, and somewhat smaller value for the H–like atoms). Even taking into account the fact that the ionized continuum includes also reflection from the warm reflector (Bianchi et al. 2001 estimated that the hot reflector contributes to about 3/4 of the total ionized reflection), from the abovementioned figures of Matt et al. (1996) a column density of the hot ionized reflector of a few$`\times 10^{21}`$ cm-2 can be deduced.

## 5 Summary

We have analysed the high energy ($`>`$4 keV) XMM–$`Newton`$ EPIC–pn spectrum of NGC 1068. The main results can be summarized as follows:

i) Possible (but not conclusive) evidence for variations of both neutral and ionized reflectors with respect to the December 1996 BeppoSAX observation is found, suggesting an upper limit to the size of both reflectors of a few parsecs. Another XMM–$`Newton`$ observation, to compare results obtained with similar instruments, is however required to confirm this finding. It is worth noting that a variability on a time scale of about a year of the ionized reflector was found by Guainazzi et al. (2000) comparing two BeppoSAX observations.

ii) Iron is overabundant, with respect to lower Z elements and when compared to the solar value, by a factor about 2. Nickel is, in its turn, overabundant by a factor $``$2 with respect to iron. A qualitatively similar result was found in the Circinus Galaxy (Molendi et al. 2003), where however both overabundances were smaller.

iii) The iron K$`\beta `$/K$`\alpha `$ flux ratio points to a low ionization state of iron, inconsistent with the K$`\alpha `$ line energy, suggesting calibration problems. The nickel K$`\beta `$/K$`\alpha `$ is larger than expected, but the K$`\beta `$ line may actually be a blend with ionized iron lines.

iv) The Fe K$`\alpha `$ Compton Shoulder is also detected, with a relative flux of 0.2, in agreement with the value expected for reflection from Compton–thick matter (Matt 2002)

v) Be–, Li–, He– and H–like iron emission lines, as well as He–like Ni line, are also found. Their EW suggests a column density of the hot ionized reflector of a few$`\times 10^{21}`$ cm-2 (Matt et al. 1996).

## Acknowledgements

Based on observations obtained with XMM–$`Newton`$, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). GM and SB acknowledge financial support from ASI.


format: text

text:
# Spin Hall effect and Berry phase in two dimensional electron gas

(October 10, 2003)

## Abstract

The spin Hall effect is investigated in a high mobility two dimensional electron system with the spin-orbital coupling of both the Rashba and the Dresselhaus types. A spin current perpendicular to the electric field is generated by either the Rashba or the Dresselhaus coupling. The spin Hall conductance is independent of the stength of the coupling, but its sign is determined by the relative ratio of the two couplings. The direction of spin current is controllable by tuning the magnitude of the surface electric field perpendicular to the two dimensional plane via adjusting the Rashba coupling. It is observed that the spin Hall conductance has a close relation to the Berry phase of conduction electrons.

Introduction.– Spintronics in semiconductor physics has become an emerging field of condensed matter because of its potential application in information industry, and also because of many essential questions on fundamental physics of electron spin.Prinz98science ; Wolf01science  It is tempting to use spin rather than charge of electron for information precessing and storage. Recently a surprising effect was predicted theoretically that an electric field can generate a dissipationless quantum spin current perpendicular to the charge current at room temperatures in conventional hole-doped semiconductors such as Si, Ge, and GaAs.Murakami03science  The effect was found to be intrinsic in electron systems with substantial spin-orbit coupling, and the spin Hall conductance has a universal value.Sinova03xxx  Based on this effect it is possible to realize spin injection in paramagnetic semiconductors rather than from the spin polarized carriers in ferromagnetic metals,Ohno99nature  which is thought as a key step to realize practical spintronic devices. The spin Hall effect due to magnetic impurities was discussed extensively.Dyakonov71 ; Hirsch99 ; Zhang00prl ; Hu03prb  When a charge current circulates in a paramagnetic metal a transverse spin imbalance will be generated, which gives rise to a spin Hall voltage or a spin current.

In two dimensional (2D) semiconductor heterostructures the spin-orbit interaction can be described in terms of two dominant contributions to the model Hamiltonian, which are directly related to symmetry of the low dimensional geometry. One type is Rashba coupling stemming from the structure inversion asymmetry of confining potential,Rashba60

$$H_R=\lambda \sigma (z\times 𝐤)=\lambda (k_x\sigma _yk_y\sigma _x)$$
(1)
where $`\sigma _a`$ ($`\alpha =x,y,z`$) are the Pauli matrices and $`k_x`$ and $`k_y`$ are two components of the wave vector. The strength of the coupling $`\lambda `$ can be modified by a gate field up to 50$`\%`$,Nitta97prl  such that it can be applied to control spin transport, such as in spin field transistor.Datta90apl ; Shen03  Another type is Dresselhaus coupling from the bulk inversion asymmetry,Dresselhaus54

$$H_D=\beta \left(k_x\sigma _xk_y\sigma _y\right).$$
(2)
In some materials such as GaAs the two types of spin-orbit coupling are usually of the same order of magnitudes. The interplay of them has been investigated theoretically with respect to several phenomena, such as the nonballistic spin field effect transistorSchliemann03prl  and electron-spin manipulationRashba03prl .

In this letter we investigate the spin Hall effect in a 2D electron system with spin-orbit coupling of both Rashba and Dresselhaus types. Consider the Hamiltonian,

$$H_0=\frac{\mathrm{}^2}{2m}𝐤^2+H_R+H_D.$$
(3)
We find the Berry phase in the eigenstates to be $`0`$ for $`\left|\lambda \right|=\left|\beta \right|`$, $`+\pi `$ for $`\left|\lambda \right|>\left|\beta \right|,`$ and $`\pi `$ for $`\left|\lambda \right|<\left|\beta \right|.`$ When the system is subjected to an electric field in the 2D plane the spin-orbit coupling leads to a dissipationless spin current perpendicular the electric field, and polarized in the direction perpendicular to the 2D plane when $`\left|\lambda \right|\left|\beta \right|.`$ The spin Hall conductivity has a universal value except that its sign is determined by the relative ratio of two couplings or the sign of the Berry phase. As Rashba coupling is tunable by a gate field perpendicular to the 2D plane it is possible to control the direction of spin current by adjusting the magnitude rather than direction of the gate field when the system is near $`\left|\lambda \right|=\left|\beta \right|.`$

Eigenstates of H0 and the Berry phase.– The system with two types of coupling has been investigated by several authors. The Hamiltonian can be diagonalized exactly. The eigenstates are

$`|k_+,+,\theta `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}e^{i\theta }\\ +i\end{array}\right)|k_+,`$ (4c)
$`|k_{},,\theta `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}e^{i\theta }\\ i\end{array}\right)|k_{}.`$ (4f)
where $`|k`$ is the eigenket of k, and $`\theta `$ is given by
$$\mathrm{tan}\theta =\frac{\lambda k_y\beta k_x}{\lambda k_x\beta k_y}=\frac{\lambda \mathrm{sin}\phi \beta \mathrm{cos}\phi }{\lambda \mathrm{cos}\phi \beta \mathrm{sin}\phi }$$
(5)
with $`k_x=k\mathrm{cos}\phi `$ and $`k_y=k\mathrm{sin}\phi .`$ We choose $`k_\pm `$ such that two eigenstates are degenerated

$$E(k)=\frac{\mathrm{}^2}{2m}k_+^2k_+\mathrm{\Delta }\omega (\phi )=\frac{\mathrm{}^2}{2m}k_{}^2+k_{}\mathrm{\Delta }\omega (\phi )$$
(6)
with $`\mathrm{\Delta }\omega \left(\phi \right)=\sqrt{\lambda ^2+\beta ^22\lambda \beta \mathrm{sin}2\phi }`$. In general the two bands do not cross over except at $`k=0`$. In the case of $`\lambda =\pm \beta ,`$ $`\theta =\pi /4`$ if $`\lambda =\beta `$ and $`\pi /4`$ if $`\lambda =\beta .`$ The spin states in Eq.(4) is independent of $`k`$ and $`\phi .`$ Furthermore $`k\mathrm{\Delta }\omega =0`$ when $`k_x=k_y`$ or $`\phi =\pi /4.`$ It is noted that there exists an additional conserved quantity, $`\left(\sigma _x\pm \sigma _y\right)/\sqrt{2}`$.Schliemann03prl  The difference of two $`k_\pm `$ is given by

$$k_+\left(\phi \right)k_{}(\phi )=\frac{2m}{\mathrm{}^2}\mathrm{\Delta }\omega (\phi )$$
(7)
which is independent of $`k.`$ The two energy bands are plotted in Fig. 1. In these states the Berry phase,Berry84  which is acquired by a state upon being transported around a loop in the k space, can be evaluated exactly,

$$\gamma _\pm =dlk_\pm ,\pm ,\theta \left|i\frac{}{𝐤}\right|k_\pm ,\pm ,\theta =\frac{\lambda ^2\beta ^2}{\left|\lambda ^2\beta ^2\right|}\pi $$
(8)
if $`\left|\lambda \right|\left|\beta \right|`$, and is equal to zero if $`\left|\lambda \right|=\left|\beta \right|.`$ The disappearance of the Berry phase at $`\left|\lambda \right|=\left|\beta \right|`$ may be relevant to the crossover of two energy bands. We notice that if we make a replacement, $`\lambda \beta `$ and $`\beta \lambda `$ the angle $`\theta \pi /2\theta .`$ This property gives different Berry phases for $`\left|\lambda \right|>\left|\beta \right|`$ and $`\left|\lambda \right|<\left|\beta \right|.`$ According to the theory of Sundaram and Niu,Sundaram99  the Berry phase is closely related to the electron transport in semiconductors. It is found that the Berry phase and magnetic monopoles in momentum space generate anomalous Hall effect in some ferromagnetic materials.Fang03

Quantum dynamics in an electric field.– We shall consider the effect of uniform electric field $`E`$. Our whole Hamiltonian is thus given by $`H=H_0+eEx`$, assuming that the electric field is along the x-direction. We start with the Heissenberg equation of motion, $`i\mathrm{}\frac{d}{dt}O=[O,H]`$. The wave vector is determined by $`k_x(t)=k_{x0}\frac{eE}{\mathrm{}}t`$ and $`k_y(t)=k_{y0}`$ where $`k_{y0}`$ and $`k_{y0}`$ are the initial values at $`t=0`$. The set of equations for spin operators are given by

$`{\displaystyle \frac{}{t}}\sigma _x`$ $`=`$ $`{\displaystyle \frac{2}{\mathrm{}}}(\lambda k_{x0}\beta k_{y0})\sigma _z+{\displaystyle \frac{2\lambda eE}{\mathrm{}^2}}t\sigma _z`$ (9a)
$`{\displaystyle \frac{}{t}}\sigma _y`$ $`=`$ $`{\displaystyle \frac{2}{\mathrm{}}}(\lambda k_{y0}\beta k_{x0})\sigma _z{\displaystyle \frac{2\beta eE}{\mathrm{}^2}}t\sigma _z`$ (9b)
$`{\displaystyle \frac{}{t}}\sigma _z`$ $`=`$ $`+{\displaystyle \frac{2}{\mathrm{}}}(\lambda k_{x0}\beta k_{y0})\sigma _x+{\displaystyle \frac{2}{\mathrm{}}}(\lambda k_{y0}\beta k_{x0})\sigma _y`$ (9c)
$`{\displaystyle \frac{2\lambda eE}{\mathrm{}^2}}t\sigma _x+{\displaystyle \frac{2\beta eE}{\mathrm{}^2}}t\sigma _y`$
When $`E=0`$, the equations can be solved exactly
$`\sigma _x^0(t)`$ $`=`$ $`\left(\mathrm{cos}^2\theta \sigma _x(0)+\mathrm{sin}\theta \mathrm{cos}\theta \sigma _y(0)\right)\mathrm{cos}\omega t`$
$`\mathrm{cos}\theta \sigma _z(0)\mathrm{sin}\omega t+\mathrm{sin}^2\theta \sigma _x(0)\mathrm{sin}\theta \mathrm{cos}\theta \sigma _y(0);`$
$`\sigma _y^0(t)`$ $`=`$ $`\left(\mathrm{sin}\theta \mathrm{cos}\theta \sigma _x(0)+\mathrm{sin}^2\theta \sigma _y(0)\right)\mathrm{cos}\omega t`$
$`\mathrm{sin}\theta \sigma _z(0)\mathrm{sin}\omega t\mathrm{sin}\theta \mathrm{cos}\theta \sigma _x(0)+\mathrm{cos}^2\theta \mathrm{cos}\theta \sigma _y(0);`$
$`\sigma _z^0(t)`$ $`=`$ $`\sigma _z(0)\mathrm{cos}\omega t+(\mathrm{cos}\theta \sigma _x(0)+\mathrm{sin}\theta \sigma _y(0))\mathrm{sin}\omega t.`$ (10c)
The characteristic frequency $`\omega =\frac{2}{\mathrm{}}k\mathrm{\Delta }\omega (\phi ).`$ Until now we have no general solution for the problem in the presence of an electric field. In this letter we are concerned with the linear response of the transport properties to the field. For our purpose we only need an asymptotic solution for a weak field $`E`$ and a shot instant $`t`$. We expand the spin operators in terms of $`E`$: $`\sigma _a(t)=\sigma _\alpha ^0(t)+\mathrm{\Delta }\sigma _\alpha ^1(t)+\mathrm{\Delta }\sigma _\alpha ^2(t)+\mathrm{}.`$ For a short instant $`t`$, we have an asymptotic solution for the linear correction to $`\sigma _\alpha ^0(t)`$ due to the electric field,
$`\mathrm{\Delta }\sigma _x^1`$ $``$ $`+{\displaystyle \frac{eE}{2k^2}}{\displaystyle \frac{\lambda \sigma _z^0(t)}{\left[\mathrm{\Delta }\omega \left(\phi \right)\right]^2}};`$ (11a)
$`\mathrm{\Delta }\sigma _y^1`$ $``$ $`{\displaystyle \frac{eE}{2k^2}}{\displaystyle \frac{\beta \sigma _z^0(t)}{\left[\mathrm{\Delta }\omega \left(\phi \right)\right]^2}};`$ (11b)
$`\mathrm{\Delta }\sigma _z^1`$ $``$ $`{\displaystyle \frac{eE}{2k^2}}{\displaystyle \frac{\lambda \sigma _x^0(t)\beta \sigma _y^0(t)}{\left[\mathrm{\Delta }\omega \left(\phi \right)\right]^2}}.`$ (11c)
We see that the effective field caused by the spin-orbit coupling has a non-zero correction to the spin operators even when $`t0.`$
Spin current and universal spin Hall conductance.– Of cause the electric field will induce a charge current in an electron gas. The problem has already been studied extensivelySchliemann03xxx . Spin-orbit coupling leads to a non-zero Hall resistance proportional to $`\lambda `$ and $`\beta ,`$ $`\rho _{xy}sign(\lambda )\left|\lambda \beta \right|.`$ In this letter we focus on the spin Hall effect. We notice that the electric field causes a linear correction to spin operators in Eqs.(11) when $`t0.`$ This correction will generate the spin Hall effect as discussed by Murakami et al.Murakami03science  and Sinova et al.Sinova03xxx . Spin current, polarized to perpendicular the 2D plane and flowing perpendicular to the electric field,Note1  is defined as

$$j_y^z=\frac{\mathrm{}}{4}\{\sigma _z(t),\frac{}{t}y(t)\}.$$
(12)
where $`\frac{}{t}y=\frac{\mathrm{}}{m}k_{y0}+\frac{\lambda }{\mathrm{}}\sigma _x(t)+\frac{\beta }{\mathrm{}}\sigma _y(t).`$ The expectation value is taken over the energy eigenstates of electrons. The linear term to the electric field is given by

$`j_y^z`$ $`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \frac{\mathrm{}^2}{2m}}k_+,+,\theta \left|\mathrm{\Delta }\sigma _z^1(t)k_{y0}\right|k_+,+,\theta n_F(E(k_+)\mu )`$ (13)
$`+{\displaystyle \underset{k}{}}{\displaystyle \frac{\mathrm{}^2}{2m}}k_{},,\theta \left|\mathrm{\Delta }\sigma _z^1(t)k_{y0}\right|k_{},,\theta n_F(E(k_+)\mu )`$
where $`n_F(E(k))`$ is the Dirac-Fermi distribution. At zero temperature we have

$$j_y^z=\frac{\mathrm{}^2eE}{16\pi ^2m}_0^{2\pi }𝑑\phi \frac{\left(\lambda ^2\beta ^2\right)\mathrm{sin}^2\phi \left[k_+k_{}\right]}{(\lambda ^2+\beta ^22\lambda \beta \mathrm{sin}2\phi )^{3/2}}.$$

By using the relation in Eq.(7) we obtain the spin Hall conductance

$$\sigma _{sH}=j_y^z/E=\frac{e}{8\pi ^2}\gamma _\pm $$
(14)
as long as both bands are occupied.Note2  We have ignored the time-dependent terms. It is an interesting observation that the spin Hall conductance is closely relevant to the Berry phase of electron states in the absence of the electric field. It deserves for further investigation on the relation in general cases. Spin current caused by the Berry phase should not contribute to dissipation, although a charge current will do. If $`\beta =0`$ $`\sigma _{sH}=e/8\pi `$ which recovers Sinova et al.’s resultSinova03xxx  and followed by several authors more recentlySchliemann03yyy ; Hu03xxx . If $`\lambda =0,`$ Dresselhaus coupling can also produce a spin Hall conductance $`\sigma _{sH}=e/8\pi `$. Very interestingly two spin Hall conductances differ from only a sign. Competition of two couplings also cancel the effect at $`\lambda =\pm \beta .`$ It seems that the effect caused by Dresselhaus coupling dominates if $`\left|\beta \right|>\left|\lambda \right|`$ while the effect caused by Rashba coupling dominates if $`\left|\lambda \right|>\left|\beta \right|`$. Very recently RashbaRashba03xxx  observed that even in the absence of an external electric field there exist a tiny spin current $`J_y^x=J_x^y0`$ in the open boundary condition. The type of spin current along spin z-direction does not exist in the absence and presence of an electric field.

To explore the physical origin of the sign change of spin Hall conductance near $`\lambda =\beta `$ we find a new symmetry in the system with the two spin-orbital couplings. Performing an unitary transformation,

$$\sigma _x\sigma _y;\text{ }\sigma _y\sigma _x;\text{ }\sigma _z\sigma _z,$$
(15)
the two terms of Rashba and Dresselhaus couplings are changed to

$`H_R`$ $`=`$ $`\lambda (k_x\sigma _yk_y\sigma _x)\lambda (k_x\sigma _xk_y\sigma _y);`$ (16a)
$`H_D`$ $`=`$ $`\beta (k_x\sigma _xk_y\sigma _y)\beta (k_x\sigma _yk_y\sigma _x).`$ (16b)
The spin current along the z-direction changes a minus sign, $`J_y^zJ_y^z`$ as $`\sigma _z\sigma _z.`$ Thus this symmetry indicates that the system with Rashba coupling $`\lambda `$ and Dresselhaus coupling $`\beta `$ and the system with Rashba coupling $`\beta `$ and Dresselhaus coupling $`\lambda `$, have the same spin Hall conductance (along spin z direction), but with opposite signs. At the symmetric point $`\lambda =\beta `$, we conclude that $`J_y^z=0`$, i.e., the spin Hall effect is suppressed completely. This general conclusion is based on symmetry analysis of the system and is mathematically rigorous. It is in agreement with the result of the linear response calculation in Eq. (14).
Now we come to address the effects of a finite electron quasi-particle lifetime and scattering mechanism, which was discussed in a system with Rashba coupling by several authors. Schliemann and LossSchliemann03yyy  discussed the correction of spin conductance due to the effect of finite lifetime $`\tau `$ for $`\beta =0`$ using the Kubo formula in the zero frequency limit and concluded that the universal value for spin conductance is valid only when the “Rashba energy” $`\epsilon _R=m\lambda ^2/\mathrm{}^2`$ should be comparable with the energy scale $`\mathrm{}/\tau .`$ Fortunately in some semiconductors such as GaAs quantum wells $`\mathrm{}/\tau <\epsilon _R`$ and it is safely to neglect these effects. On the other hand, more recently Burkov and MacDonald found that the spin Hall conductance in a system with Rashba coupling remains the universal value $`e/8\pi `$ in the case of strong disordering as in the case of weak disordering.Burkov03xxx  Thus the scattering effect depends on the scattering mechanisms.

Controllable Rashba coupling and revertible spin current.– Rashba coupling and Dresselhaus coupling are related to the symmetry of quantum wells. In GaAs quantum wells both terms are usually of the same order of magnitude while the narrow gap compounds like in InAs Rashba coupling dominates.Jusserand95  However there is no difficulty to achieve the situation of $`\lambda =\beta `$. In the case of $`\lambda =\pm \beta `$ the Elliot-Yafet spin flip mechanism is suppressed completely since the spin states in $`|k_\pm ,\pm ,\theta =\pm \pi /4`$ is independent of $`\phi `$. The Berry phase of the states $`|k_\pm ,\pm ,\theta `$ becomes zero and correspondingly there is no spin Hall effect. Thus the electric field cannot generate a dissipationless spin current. If $`\left|\lambda \right|>\left|\beta \right|`$ $`\sigma _{sH}=e/8\pi `$ while $`\sigma _{sH}=e/8\pi `$ if $`\left|\lambda \right|<\left|\beta \right|.`$ The spin current changes a sign (or direction). As Rashba coupling is adjustable by a gate field perpendicular to the electron gas planeNitta97prl  we may change $`\gamma `$ around $`\beta ,`$ $`\gamma =\beta +\delta \mathrm{\Delta }E_{},`$ by adjusting the perpendicular electric field $`\mathrm{\Delta }E_{}`$. Thus by using the property of sign change in the Berry phase we can realize to revert the spin current by adjusting the magnitude rather than the direction of the gate field.

Conclusion.– The spin Hall conductance caused by Dresselhaus coupling has the same universal value, but an opposite sign to that caused by Rashba coupling. It is an interesting observation that the spin Hall conductance has a close relation to the Berry phase of single electron states. The spin current in this effect can be revertible when the amplitudes of two couplings are comparable by adjusting the gate voltage via Rashba coupling. It is anticipated that this effect will be applicable to the spintronic devices in the future.

The author thanks S. Murakami, Q. Niu, J. Sinova, X. C. Xie, and F. C. Zhang for discussions. This work was supported by a grant from the Research Grant Council of Hong Kong (Project No.: HKU 7088/01P). After this paper was submitted for publication the author was aware that Sinova et al. arrived at the same spin Hall conductances by means of the linear response theory independently.Sinitsyn03xxx


text:
# The apsidal corotation in mean motion resonance: the 55 Cancri as an example

## 1 Introduction

Over 100 extra-solar planets have been found (see e.g. http://www.obspm.fr/encycl/encycl.html by J.Schneider). The high eccentricities, close distances to the central stars and the heavy masses, show the wide variety of planetary system that are quite different from our solar system.

Among these “exoplanets”, 28 are located around the main sequence stars in 13 “multiple planet systems”, and some of them appear to be in Mean Motion Resonances (hereafter MMR). Recently, the two inner planets around the 55 Cancri were found to be in a possible 3:1 MMR \[Ji et al. 2003\]. In numerical simulations, the resonant angles are found to librate around certain values. At the same time, the difference between the two periastrons is found to be locked to a fixed value, that is, the two planets precess at the same rate.

It is commonly accepted that these extra-solar planets are not formed in situ but have suffered orbital migrations. Different migration speeds of planets may lead to a variation of mutual distance between two neighboring planets and eventually a capture into the current resonant configurations \[Kley 2000, Kley 2003, Nelson & Papaloizou 2002\]. The dynamics of these resonances may contain important information of the capturing process and therefore give valuable hints to the planet formation and the evolution of a planetary disk. As a result, planetary systems with resonance require special attention. Moreover, the possibility of an Earth-type (habitable) planet in the 55 Cancri system has been discussed \[Cuntz et al. 2003\]. The existence of a habitable planet in a system depends not only on the current dynamical features of the system but also on the history it has experienced. Therefore it is important to study the now-known planets before an Earth-like planet was observed someday. All these points make this planetary system interesting.

In this paper, we confirm the possibility that the two planets are in a 3:1 MMR by numerical integrations and discuss the possible configurations this system would take (section 2), then with a new expansion of the Hamiltonian for a planar three-body problem (section 3), we analyse the dynamics of this system. We will discuss the occurrence of the apsidal corotation and its influence on the stability of this system (section 4). Finally the conclusions and discussions are given in section 5.

## 2 Numerical simulations

### 2.1 Numerical model and initial set-up

The 55 Cancri system is a system with three planets. We adopt the orbital solution given by Fischer et al. \[Fischer et al. 2003\], as being listed in Table 1. To numerically simulate this system, we use a symplectic integrator \[Wisdom & Holman 1991, Mikkola & Palmer 2000\], which allows us to follow the orbital evolution of each body and simultaneously the stability of this orbit indicated by the Lyapunov Character Indicator (LCI) \[Froeschlé 1984\]. Each planet orbit has an LCI and we choose the largest one among the three to indicate the stability of this planetary system in this paper. The time step is set to be $`0.3`$ d, which is about $`2\%`$ of the orbital period of the innermost planet.

We adopt the masses, semimajor axes and eccentricities of planets the values listed in Table 1. The initial orbital inclination of companion D is set to be zero, while the companion B and C have initial inclinations of $`10^5`$ degrees. Other angles (the ascending node, the periastron and the mean anomaly) are randomly generated from $`[0,2\pi )`$. Four hundred simulations, with different initial conditions, are integrated up to $`10^6`$ years.

Ignoring the companion D, we first integrate a system composed of the central star and the companion B and C. Then, retaining the initial conditions of these three bodies and adding the companion D on an initial orbit with the given $`a,e,i`$ and randomly selected angle elements, we integrate the four-body system. As the simulations show, the companion D has nearly no influence on the motion of the inner planets, although it is much heavier. This is due to its large semimajor axis ($`5.461`$ au). However, in order to get reliable simulations for the real system, we report here the numerical results from the four-body model.

During the integrations, if the distance between any two of the planets became smaller than half of the criterion for the “Hill stability” \[Gladman 1993\], the system was considered become collapsed and the simulation was terminated.

Generally, the motion of a planet is dominated by the central star and its orbit is a conic section with small deviations due to the gravitational perturbation from other planet(s). This perturbation can be described by the disturbing function, which can be expanded in terms of the orbital elements. We use $`a,e,i,\varpi ,\mathrm{\Omega },\lambda `$ to denote the semi-major axis, eccentricity, inclination, longitude of periastron, longitude of ascending node, and mean longitude, respectively. The disturbing function for a planet with mass $`m_1`$ and orbital elements $`(a_1,e_1,i_1,\varpi _1,\mathrm{\Omega }_1,\lambda _1)`$, perturbed by another planet (indicated by subscripts ‘2’), can be written as

$$R=\mu _2F(a_1,a_2,e_1,e_2,i_1,i_2)\mathrm{cos}\phi .$$

Here $`\mu _2=𝒢m_2`$ and $`\phi `$ is a linear combination of angles

$$\phi =j_1\lambda _1+j_2\lambda _2+j_3\varpi _1+j_4\varpi _2+j_5\mathrm{\Omega }_1+j_6\mathrm{\Omega }_2$$

where $`j_1,j_2,\mathrm{},j_6`$ are integers satisfying $`j_1+j_2+\mathrm{}+j_6=0`$.

Particularly, when the orbital periods satisfy $`P_2/P_13`$, that is, $`\dot{\lambda }_13\dot{\lambda }_20`$ (dots denote the time derivation), we have $`\dot{\phi }_{(j_1=1,j_2=3)}0`$ since generally $`\dot{\varpi }_1,\dot{\varpi }_2,\dot{\mathrm{\Omega }}_1,\dot{\mathrm{\Omega }}_2\dot{\lambda }_1,\dot{\lambda }_2`$. Then in the disturbing function $`R`$, those terms containing $`\lambda _13\lambda _2`$ can not be eliminated by the averaging technique and become the leading terms. And if the inclinations $`i_1,i_2`$ are not too large, the leading terms in $`R`$ are of $`O(e_1^2),O(e_2^2),O(e_1e_2)`$ and the corresponding angle $`\phi `$, now called resonant arguments, are $`\theta _1=\lambda _13\lambda _2+2\varpi _1`$, $`\theta _2=\lambda _13\lambda _2+2\varpi _2`$, $`\theta _3=\lambda _13\lambda _2+\varpi _1+\varpi _2`$. Therefore the libration of any one of $`\theta _{1,2,3}`$ indicates the two planets are in an (eccentricity-type) 3:1 MMR. We note that no more than two of them are linear independent and $`\theta _3=(\theta _1+\theta _2)/2`$ \[see for example Murray & Dermott \[Murray & Dermott 1999\] for more details\].

Beside these resonant arguments, the relative apsidal longitude $`\mathrm{\Delta }\varpi =\varpi _1\varpi _2`$ is also a critical argument often to be discussed. We call the libration of $`\mathrm{\Delta }\varpi `$ in an MMR the apsidal corotation. In a 3:1 MMR, if both $`\theta _1`$ and $`\theta _2`$ librate with small amplitudes, $`\mathrm{\Delta }\varpi =\frac{1}{2}(\theta _1\theta _2)`$ will also necessarily librate. In this sense, the apsidal corotation is not independent and should be different from the term of apsidal resonance, which is in the context of a secular perturbation.

Hereafter we use subscripts $`1,2`$ to label the orbital elements of companion B and C respectively.

### 2.2 Numerical results

About one third ($`133`$ out of $`400`$) of the simulations collapse during the integrations, and all the remainders survive for the $`10^6`$ yr integration. Most of the survivors have very short ($`<10^2`$ yr) e-folding time $`T_e`$ ($`=1/\mathrm{LCI}`$, the reciprocal of LCI), while 38 of them have $`T_e10^3`$ yr, and they are in this sense regarded as stable. In these stable systems, all the three planets have final inclinations $`0\mathrm{°}.1`$, but those unstable ($`T_e<10^3`$ yr) systems may have inclinations of companion B and C as high as $`20\mathrm{°}`$, that is, the stable system prefers to retain to be coplanar. All the stable systems are found to be associated with the 3:1 MMR between the two inner planets.

According to different configurations of $`\theta _{1,2,3}`$ and $`\mathrm{\Delta }\varpi `$, the stable systems can be divided into three groups, with representative cases illustrated in Fig. 1. In all the cases, the companion D is in a steady motion (we don’t illustrate this here) and, the semi-major axes of companion B and C vibrate with very small variations (an example is shown in Fig. 2).

Case a (Fig. 1) shows a strong resonance with $`\theta _1`$, $`\theta _2`$ and $`\theta _3`$ librating around $`215\mathrm{°}`$, $`75\mathrm{°}`$ and $`325\mathrm{°}`$, respectively. Simultaneously, the $`\mathrm{\Delta }\varpi `$ also librates around $`250\mathrm{°}`$. The well known symmetric configurations of $`\mathrm{\Delta }\varpi `$ librating around $`0\mathrm{°}`$ (alignment) or $`180\mathrm{°}`$ (anti-alignment) imply the conjunctions of the two planets always occur at a certain position and this attributes to the stability of the system \[Lee & Peale 2002\]. However the asymmetric configuration with $`\mathrm{\Delta }\varpi `$ librating around angles other than $`0\mathrm{°}`$ or $`180\mathrm{°}`$, has been predicted by Beaugé, Ferraz-Mello & Michtchenko \[Beaugé et al 2003\].

Case b differs from case a in the behavior of $`\theta _1`$, $`\theta _3`$ and $`\mathrm{\Delta }\varpi `$, which circulate now. The $`\theta _2`$ still librates around $`80\mathrm{°}`$ with an amplitude of $`40\mathrm{°}`$. In fact, when $`\theta _2`$ librates with a small amplitude, we may assume $`\theta _2const.`$, that is, $`\dot{\theta }_20`$. From $`\dot{\theta }_2=\dot{\lambda }_13\dot{\lambda }_2+2\dot{\varpi }_20`$ we can derive

$$\frac{n_2\dot{\varpi }_2}{n_1\dot{\varpi }_2}\frac{1}{3},$$

where $`n_{1,2}=\dot{\lambda }_{1,2}`$ are the mean motions of planets. In a reference frame corotating with the periastron of the second planet, the relative mean motions are $`n_k^{}=n_k\dot{\varpi }_2`$ $`(k=1,2)`$, and from $`n_2^{}/n_1^{}1/3`$ we know that in this reference frame the conjunctions of planets always happen at a certain direction. Meanwhile, given $`n_k\dot{\varpi }_2`$, we have $`n_2/n_1n_2^{}/n_1^{}`$ and there is still a 3:1 relation between the two orbital periods. In the 55 Cancri system, our numerical results show $`n_{1,2}`$ are larger than $`\dot{\varpi }_2`$ by at least two order of magnitudes.

Although the numerical simulations give more examples with configuration as case a than case b (28 vs 7, see Table 2), we notice that among the most stable sets with e-folding time $`T_e10^6`$ yr (the same order as the integration time), there are two sets having configurations as case b but only one as case a. This prevents us from concluding which type, case a or case b, is more stable.

In case c, $`\theta _1,\theta _3`$ and $`\mathrm{\Delta }\varpi `$ circulate, while $`\theta _2`$ librates around $`180\mathrm{°}`$ with a quite large amplitude of $`160\mathrm{°}`$. As a result, it is less stable than the above two cases. The only three examples with this configuration in our simulations all have e-folding times $`<5\times 10^3`$ yr.

Finally we summarize the e-folding times and the motion configurations of the 38 stable systems in Table 2.

Besides the angles $`\theta _{1,2,3}`$ and $`\varpi `$, we notice that the orbital eccentricities $`e_1,e_2`$ also have remarkably different behaviors in these three cases. Another noteworthy feature of numerical results is that the system has a symmetry over angles $`\theta _i`$ and $`\mathrm{\Delta }\varpi `$, that is, when $`\varphi `$ is a libration center of $`\theta _i`$ or $`\mathrm{\Delta }\varpi `$, the argument $`(360\mathrm{°}\varphi )`$ is too. For example, we observe $`\theta _2`$ librates around $`80\mathrm{°}`$ and $`280\mathrm{°}`$ with nearly the same ratio in our simulations.

### 2.3 Some extended calculations

To integrate all the surviving systems to $`10^7`$ yr or even longer is out of our current computational capability. However, to explore further the stabilities of different configurations, we integrate the three examples shown in Fig. 1 further to $`3\times 10^7`$ years. Both the case a and b last this time span and keep the configurations all the time. But case c loses the stability when the inclinations of companion B and C increase to $`5\mathrm{°}`$ after $`1\times 10^7`$ years. As in those unstable examples, this inclination increasing is due to a sequence of close encounters between the two planets. On the contrary, the mechanism of the 3:1 MMR can protect those systems as case a and b from catastrophic close encounters.

Other 5 examples randomly selected from the 17 systems with $`T_e10^3`$ yr are integrated to $`1\times 10^7`$ yr. Again, except the one having motion configuration as case c, they all survive and have new e-folding time $`T_e>10^3`$ yr. \[The LCI will converge to a certain value as the integrating time approaches to infinity, so $`T_e`$ of an orbit would change when it’s integrated to a longer time span. see for example \[Froeschlé 1984\].\] As for those systems with larger $`T_e`$, they are more stable and should have longer surviving lives.

So far, our discussions are all based on the small initial inclinations of the two inner planets. We find that in stable systems the inner two planets are locked to the 3:1 MMR and their small orbital inclinations can be kept to at least $`10^6`$ yr. While in those unstable systems, chaotic close encounters of planets could cause inclinations to increase. Furthermore, we also numerically simulated hundreds of examples with larger initial inclinations. In these systems, except the inclinations of the two inner planets are set to higher values, all the orbital elements are initialized in the same way as before, and we integrate these systems to $`10^6`$ yr. A brief summary is given in Table 3.

We see the survive probability drops from $`66\%`$ to 21% as the initial inclinations increase from $`10^5`$ to 25 degrees, but at the same time the percentage of stable systems even does not decrease until $`i_{1,2}^0=10\mathrm{°}`$ before it attains a value of $`5\%`$. This implies that the stability of the system is not so sensitive to the inclinations, and if smaller than 10 degrees, they have nearly no influence on the stability of the system. On the other hand, that all the stable examples are associated with the 3:1 MMR reveals again the importance of this resonance in stabilizing this planetary system, even when the orbits are far from coplanar. We don’t give more analysis to the problem of higher inclinations here, but it deserves a detailed discussion in future .

## 3 Analytical model

To understand the dynamics of the 3:1 MMR in the 55 Cancri, it is useful to describe it with an analytical model. As we mentioned above, the companion D has a very weak influence on the motion of the inner planets. In addition, the numerical simulations show it is mainly the 3:1 MMR that is responsible for the stability of the system, in spite of the higher orbital inclinations. So, a planar three-body problem is a reasonable model for describing this planetary system, especially when we mainly focus on the 3:1 MMR. Of course, considering the possible substantial difference between the coplanar and non-coplanar motion, a more realistic and reliable model is still needed. On the other hand, the usual analytical perturbation models, such as the Laplace expansions, are not convergent for high eccentricities, and expansions with low-order truncations may yield imprecise results \[Ferraz-Mello 1994\]. Recently, a new expansion of the Hamiltonian suitable for high-eccentricities has been developed by Beaugé and Michtchenko \[Beaugé & Michtchenko 2003\]. We will apply this expansion to the system the 55 Cancri.

### 3.1 The Hamiltonian

Suppose three bodies with mass $`M_0,m_1,m_2`$ orbiting around their common center of mass. $`M_0m_{1,2}`$ is the central star and $`m_{1,2}>0`$ are the planets. Near the $`3:1`$ MMR, a set of canonical variables are defined as:

$`\lambda _1;J_1=L_1+{\displaystyle \frac{1}{2}}(I_1+I_2)`$
$`\lambda _2;J_2=L_2{\displaystyle \frac{3}{2}}(I_1+I_2)`$
$`\sigma _1={\displaystyle \frac{3}{2}}\lambda _2{\displaystyle \frac{1}{2}}\lambda _1\varpi _1;I_1=L_1G_1`$ (1)
$`\sigma _2={\displaystyle \frac{3}{2}}\lambda _2{\displaystyle \frac{1}{2}}\lambda _1\varpi _2;I_2=L_2G_2`$
where $`L_i=m_i^{}\sqrt{\mu _ia_i}`$ and $`G_i=L_i\sqrt{1e_i^2}`$ are Delaunay variables, $`a_i,e_i,\lambda _i,\varpi _i`$ are the orbital elements of the $`i`$th planet, $`\mu _i=𝒢(m_i+M_0)`$, $`𝒢`$ is the gravitational constant, and $`m_i^{}=\frac{m_iM_0}{m_i+M_0}`$ is the reduced mass. These canonical variables define a four degree-of-freedom system. And let’s bear in mind $`\theta _i=2\sigma _i`$ ($`i=1,2`$), $`\theta _3=\sigma _1\sigma _2`$ and $`\mathrm{\Delta }\varpi =\sigma _2\sigma _1`$.

After considering the conservation of angular momentum and averaging over the fast-angle, we get two integrals of motion $`J_{\mathrm{sum}}=J_1+J_2`$, $`J_{\mathrm{dif}}=J_1J_2`$. At the same time, the problem is reduced to a two degree-of-freedom system in canonical variables ($`I_i,\sigma _i`$). The resulting Hamiltonian reads

$$H=\underset{i=1}{\overset{2}{}}\frac{\mu _i^2m_i^3}{2L_i^2}+H_1$$
(2)
where $`H_1`$ is the truncated disturbing function,

$`H_1={\displaystyle \frac{𝒢m_1m_2}{a_2}}{\displaystyle \underset{n=0}{\overset{3}{}}}{\displaystyle \underset{j=0}{\overset{j_{\mathrm{max}}}{}}}{\displaystyle \underset{k=0}{\overset{k_{\mathrm{max}}}{}}}{\displaystyle \underset{u=0}{\overset{u_{\mathrm{max}}}{}}}{\displaystyle \underset{l=l_{\mathrm{max}}}{\overset{l_{\mathrm{max}}}{}}}R_{n,j,k,u,l}(\alpha \alpha _0)^n`$
$`\times e_1^je_2^k\mathrm{cos}(2u\sigma _1+l(\sigma _2\sigma _1)).`$ (3)
Here $`\alpha =a_1/a_2`$ is the ratio between the two semimajor axes and $`\alpha _0`$ is its value in the exact MMR. The coefficients $`R_{n,j,k,u,l}`$ can be determined beforehand \[Beaugé & Michtchenko 2003\]. Considering the high value of eccentricities, we adopt the upper limits of sums $`u_{\mathrm{max}}=l_{\mathrm{max}}=12`$ and $`j_{\mathrm{max}}=k_{\mathrm{max}}=15`$ in this paper. For convenience, the Hamiltonian has been expressed explicitly in $`(a_i,e_i)`$ instead of $`(I_i,J_i)`$, therefore the integrals $`J_{\mathrm{sum}}`$ and $`J_{\mathrm{dif}}`$ don’t appear explicitly but exist as constrains. We note the total energy $`H`$ is an integral too. Furthermore, it is easy to prove the symmetry buried in the Hamiltonian:

$$H(\sigma _1,\sigma _2)H(\sigma _1,\sigma _2)H(180\mathrm{°}+\sigma _1,180\mathrm{°}+\sigma _2).$$

### 3.2 The reliability of the Hamiltonian

We can calculate planet orbits by applying this Hamiltonian. Since it is an averaged system, the time step can now be much longer than that used in our symplectic integrator. A comparison between this Hamiltonian expansion and the direct numerical integration is shown in Fig. 3, in which we choose the most unstable case of Fig. 1c as an example. We get the same initial conditions and recalculate the orbits from the Hamiltonian equations with an average (variable) time step of $`60`$ d, which is four times of the inner planet’s orbital period.

The agreements in Fig. 3, not only imply the reliability and suitability of this Hamiltonian, but also reflect that the third planet can be ignored when studying the dynamics of the inner two planets.

## 4 The apsidal corotation and the stability of the system

The Hamiltonian introduced above has two degrees of freedom. Although it is still non-integrable, it can help us to understand the dynamical properties of the system. In this section, we will use it to find out why and when the system follows the different types of resonant evolutions as shown in Fig. 1, and then via the surface of section technique we study the stability of the systems with different configurations.

### 4.1 The eccentricities and the resonant angles

Since $`J_{\mathrm{sum}}`$ and $`J_{\mathrm{dif}}`$ are integrals, $`J_1=\frac{1}{2}(J_{\mathrm{sum}}+J_{\mathrm{dif}})`$ and $`J_2=\frac{1}{2}(J_{\mathrm{sum}}J_{\mathrm{dif}})`$ are invariables too. They constrain the variation of the planetary semi-major axes by $`3J_1+J_2=3m_1^{}\sqrt{\mu _1a_1}+m_2^{}\sqrt{\mu _2a_2}=const.`$ However, when we discuss planets trapped in an MMR, we may simply assume the semi-major axes are constants, otherwise the resonance would be destroyed. In fact, the numerical simulations in section 2 also support this assumption (Fig. 2). Similarly, simple calculations show that

$$J_{\mathrm{sum}}=m_1^{}\sqrt{\mu _1a_1(1e_1^2)}+m_2^{}\sqrt{\mu _2a_2(1e_2^2)}.$$
(4)
This confines the variation of the eccentricities when we fix the semi-major axes $`a_1`$ and $`a_2`$. Besides this, the variation range of the eccentricities is also bounded by other constrains, e.g., the total energy of the system.

We assume the two semi-major axes are fixed, $`a_1=0.115`$ au and $`a_2=0.241`$ au. Then Eq. (4) defines a family of curves on the $`(e_1,e_2)`$-plane (Fig. 4). Each curve in Fig. 4 corresponds to a definite value of the integral $`J_{\mathrm{sum}}`$. The orbital elements in Table 1 give $`J_{\mathrm{sum}}=3.591733\times 10^4`$, so that the eccentricities $`e_1,e_2`$ of orbits beginning from these initial conditions are restricted on the corresponding curve during the dynamical evolution.

After assuming $`a_{1,2}`$ fixed, for certain $`e_{1,2}`$ values, we can calculate the energy level curves of the Hamiltonian on the $`(\sigma _1,\sigma _2)`$-plane. Fig. 5 shows such contours. At the first glance, we note the symmetry over both $`\sigma _1`$ and $`\sigma _2`$. This has been observed in our numerical simulations (recalling $`\theta _i=2\sigma _i`$).

Since the total energy $`H`$ is an integral, a system starting from a point on any curve in Fig. 5 can never jump to other curve with a different energy. On the other hand, the variation of eccentricities (Fig. 4) is constrained by the integral $`J_{\mathrm{sum}}`$ through Eq. (4). Therefore, as the eccentricities vary following Eq. (4), an energy curve on the $`(\sigma _1,\sigma _2)`$-plane will extend to an energy surface (e.g. a cylinder) and the system will evolve on this surface. To illustrate this, we fix $`a_1=0.115`$ au and $`a_2=0.241`$ au, set $`e_1=0.03`$ (and sequently $`0.10,0.15,0.20`$), evaluate the corresponding $`e_2`$ from Eq. (4) with $`J_{\mathrm{sum}}=3.591733\times 10^4`$, then calculate contours for different pairs of $`(e_1,e_2)`$. Finally, we get a series of section of the energy surface in Fig. 5. In this way, we can follow the dynamical evolution of a system starting from definite initial conditions.

For example, the thick curves in Fig. 5, representing an energy $`H=3.957679\times 10^3`$, keep the closed oblate shape when $`e_1`$ evolves from the initial value $`0.03`$ to $`0.20`$. The oblate profile of the curve implies that the angle $`\sigma _1`$ ($`\theta _1`$) has a larger reachable range than $`\sigma _2`$ ($`\theta _2`$)during the temporal evolution of a system with initial conditions located on this curve. The angle $`\sigma _1`$ could vary in a range of $`(160\mathrm{°},400\mathrm{°})`$, while $`\sigma _2`$ is bounded between $`20\mathrm{°}`$ and $`60\mathrm{°}`$. Correspondingly, the resonant angle $`\theta _1`$ and $`\theta _2`$ could vary in ranges with widths of $`480\mathrm{°}`$ (circulation) and $`80\mathrm{°}`$ (libration), respectively. This is just the situation of case b in Fig. 1. Another example indicated by a dashed curve with $`H=3.957722\times 10^3`$ exhibits the possibility of circulating of both $`\theta _1`$ and $`\theta _2`$.

We denote the possible variation of $`\sigma _i`$ in Fig. 5 by $`\mathrm{\Delta }\sigma _i`$. According to the foregoing analysis, $`\mathrm{\Delta }\sigma _2<\mathrm{\Delta }\sigma _1`$. If $`\mathrm{\Delta }\sigma _1180\mathrm{°}`$, both $`\theta _1`$ and $`\theta _2`$ would librate ($`\mathrm{\Delta }\theta _2<\mathrm{\Delta }\theta _1360\mathrm{°}`$) and the system could run as Fig. 1a. It could run as Fig. 1b, if $`\mathrm{\Delta }\sigma _1180\mathrm{°}`$ while $`\sigma _2`$ is still bounded in a narrow range, say $`\mathrm{\Delta }\sigma _290\mathrm{°}`$. And if $`90\mathrm{°}\mathrm{\Delta }\sigma _2180\mathrm{°}`$, Fig. 1c could be followed. So the measures of such areas in the whole $`(\sigma _1,\sigma _2)`$-plane tell the probabilities of different types of motion. Calculations show the proportions of these regions in the whole plane, i.e. the probabilities of motion type a, b and c, approximatively are $`10\%,16\%`$ and $`4.5\%`$, respectively. They are coarse estimates because we don’t count the stability of the a possible motion. After considering the stability, these estimates must decrease in different extents.

In a word, given $`a_i,e_i`$, the dynamical evolution of a system is determined by the initial values of $`\sigma _i`$ ($`i=1,2`$), which consequently correspond to different energy $`H`$.

### 4.2 The occurrence of apsidal corotation

When two planets are in a 3:1 MMR, at least one of the three resonant angles librates around a definite value. Since $`\sigma _2`$ generally has a smaller libration amplitude than $`\sigma _1`$, we can further simplify the Hamiltonian by safely setting $`\sigma _2\sigma _2^0`$. To determine $`\sigma _2^0`$, we can calculate the stable stationary solutions of the Hamiltonian equations \[Beaugé et al 2003\], or we can get an estimate from the position of the Hamiltonian extremum in Fig. 5. But in this paper, we simply adopt $`\sigma _2^0=320\mathrm{°}`$ ($`\theta _2=80\mathrm{°}`$) from the numerical results. Recalling the constant semimajor axes and the constrain on eccentricities, the Hamiltonian can be explicitly expressed as $`H=H(a_2^0,a_1^0,e_2,e_1(J_{\mathrm{sum}},e_2),\sigma _2^0,\mathrm{\Delta }\varpi )`$, with only two “free” variables $`e_2`$ and $`\mathrm{\Delta }\varpi `$. We calculate the contour of such Hamiltonian on the $`(e_2,\mathrm{\Delta }\varpi )`$-plane and show it in Fig. 6. The contour of the velocity of $`\mathrm{\Delta }\varpi `$ is also calculated and shown in Fig. 7.

A system approximately satisfying the assumptions mentioned above will evolve along one of the curves in Fig. 6. For the inner two planets in the 55 Cancri, the initial orbital elements listed in Table 1 ($`e_1^0=0.03,e_2^0=0.41`$) correspond to the right boundary of the box in Fig. 6. Two different fates of a system starting from these initial conditions are distinguished: one is characterized by a libration of $`\mathrm{\Delta }\varpi `$ around $`250\mathrm{°}`$ and the other a circulation of $`\mathrm{\Delta }\varpi `$. In the latter case, the eccentricity $`e_2`$ can reach a smaller value (correspondingly $`e_1`$ a higher value). These agree pretty well with the directly numerical integrations. Such agreements are shown in Fig. 6 by a comparing plots of the numerical results.

Fig. 7 shows the velocity of $`\mathrm{\Delta }\varpi `$. Through this figure, we can explain how $`\mathrm{\Delta }\varpi `$ evolves. Along the right boundary of the box, the absolute value of velocity could be very large, so that $`\mathrm{\Delta }\varpi `$ increases or decreases quickly. Let us write the perturbing part of the Hamiltonian as $`H_1=_{j=0}^{j_{\mathrm{max}}}_{k=0}^{k_{\mathrm{max}}}R_{j,k}e_1^je_2^k`$, where $`R_{j,k}`$ contains the coefficients $`R_{n,j,k,u,l}`$ and sums over $`n,u,l`$ in Eq.(3). By careful calculations, we obtain

$`{\displaystyle \frac{d\mathrm{\Delta }\varpi }{dt}}={\displaystyle \frac{\sqrt{1e_1^2}}{e_1^2L_1}}{\displaystyle \underset{j=1}{\overset{j_{\mathrm{max}}}{}}}{\displaystyle \underset{k=0}{\overset{k_{\mathrm{max}}}{}}}jR_{j,k}e_1^je_2^k`$
$`{\displaystyle \frac{1e_2^2}{e_2^2L_2}}{\displaystyle \underset{j=0}{\overset{j_{\mathrm{max}}}{}}}{\displaystyle \underset{k=1}{\overset{k_{\mathrm{max}}}{}}}kR_{j,k}e_1^je_2^k.`$ (5)
Obviously, when $`e_i1`$, the dominating term on the right hand side has the order $`\frac{1}{e_i}`$. Near the right boundary of the box in Fig. 7, we have $`e_10.031`$. That’s the reason why $`\mathrm{\Delta }\varpi `$ has a quick change when $`e_1`$ approaches its minimum value as shown in Fig. 1b,c.

The numerical simulations reveal that the quick varying of $`\mathrm{\Delta }\varpi `$ in Fig. 1b,c is caused by a quick varying of $`\varpi _1`$ when $`e_10`$. In fact, when an orbit is nearly circular ($`e0`$), it would be relatively easy to change its direction ($`\varpi `$), thus the two planets could loose the locking of periastrons if they were not in a strong coupling. On the other hand, whether they can strongly couple with each other depends also on how close they can approach each other, that is, how large their eccentricities (especially $`e_2`$) are. Bearing this in mind we can understand why in case a $`e_2`$ has a higher lower-limit than in cases b and c. Of course, a small $`e_2`$ can also cause a quick variation of $`\mathrm{\Delta }\varpi `$, but this situation does not happen because the energy integral prevents $`e_2`$ from approaching very small value.

From Fig. 6 and Fig. 7, we see different behaviors of $`\mathrm{\Delta }\varpi `$ (circulation or libration) have different varying directions along the right edge of the box. By numerically solving the equation $`d\mathrm{\Delta }\varpi /dt=0`$, we find two zero velocity points on the boundary. They define the boundaries of the initial conditions leading to a librating or circulating of $`\mathrm{\Delta }\varpi `$. The calculations show that the apsidal corotation happens if $`\mathrm{\Delta }\varpi ^0(160\mathrm{°},330\mathrm{°})`$, while if $`\mathrm{\Delta }\varpi ^0(30\mathrm{°},160\mathrm{°})`$, $`\mathrm{\Delta }\varpi `$ circulates.

The case c has a Hamiltonian value a little further from the extremum in Fig. 5, so that $`\sigma _2(\theta _2)`$ has a large amplitude of libration and can not be assumed as a constant. As a result, we can not find a corresponding curve for it in Fig. 6. However, if we draw a series of Hamiltonian contour with different value of $`\sigma _2`$ (adding another dimension to Fig. 6), we can also understand its dynamical evolution.

Fig. 8 summarizes the initial conditions of the 38 stable examples in the numerical simulations. For clarity, we have mapped the symmetric configurations to one regime according to the symmetry of the Hamiltonian. Except the three examples as case c, all the points for stable systems gather along the line of $`\sigma _2^0=320\mathrm{°}`$, and points for librating and circulating $`\mathrm{\Delta }\varpi `$ (case a and b) are separated by the dotted lines of $`\mathrm{\Delta }\varpi ^0=160\mathrm{°},330\mathrm{°}`$. These agree very well with the above analysis with the Hamiltonian.

We have given an estimate of the probability of following different ways as case a, b and c. The analyses here tell that all the initial conditions leading to stable configurations gather in a narrow strip along the line in Fig. 8. That is, for systems having given $`a_{1,2},e_{1,2}`$ as in Table 1, only those with initial conditions satisfying $`\sigma _2^0320\mathrm{°}`$ and $`\sigma _1^0+\mathrm{\Delta }\varpi ^0320\mathrm{°}`$ can survive for long time.

### 4.3 Surfaces of section

We know now the initial conditions determine whether the apsidal corotation happens or not, but whether a system would follow a definite evolution type also depends on the stability of an evolution type. The stability of an orbit in the phase space can be revealed by its projection on the surface of section. We present in Fig. 9 the sections of $`(I_2\mathrm{cos}\sigma _2,I_2\mathrm{sin}\sigma _2)`$ for the three orbits in Fig. 1. The plane is defined as $`\sigma _1=270\mathrm{°},\dot{\sigma }_1>0`$. To obtain these, we adopt the initial conditions of orbits from Fig. 1, and then integrate the corresponding Hamiltonian equations based on the expansion as Eq. (2,3).

Fig. 9c indicates case c has a chaotic orbit. On the contrary, the invariant curves for cases a and b (Fig. 9a,b) imply that both these two orbits are regular. That is, no matter whether the apsidal corotation happens or not, the system trapped in the 3:1 MMR could be stable. From this point of view we may argue the apsidal corotation only has a limited contribution to the stability of the system.

Because the Hamiltonian is expressed explicitly in $`a_i,e_i`$ and the integrals of motion $`J_{\mathrm{sum}}`$, $`J_{\mathrm{dif}}`$ affect the system as constrains, we can only define an energy surface ($`HH_0`$) with certain pre-determined $`J_{\mathrm{sum}},J_{\mathrm{dif}}`$, then calculate the surface of section. We find the structure of the phase space and the stability of the system depends sensitively not only on $`H`$ but also on $`J_{\mathrm{sum}}`$ and $`J_{\mathrm{dif}}`$. For example, with the same values of $`H,J_{\mathrm{sum}}`$ and $`J_{\mathrm{dif}}`$ as the ones in Fig. 1b, we calculate tens of orbits with different initial conditions. On the surface of section, these orbits occupy a big regular-motion region. This result does not conflict with the fact that there are only a small fraction of stable systems in numerical simulations, in fact, when we change $`H`$, (or $`J_{\mathrm{sum}}`$, $`J_{\mathrm{dif}}`$) just a little, we find remarkable chaotic region on the surface of section. The application of the surface-of-section technique is limited by the high dimension of the system. This will be discussed in detail in our future paper.

## 5 Summary and Discussion

With hundreds of numerical simulations of the planetary system of the 55 Cancri, we find the third planet has a very weak influence on the motion of the inner two planets. We confirm the inner two could be trapped in a 3:1 mean motion resonance and three different types of motion are found. Judging from the Lyapunov character indicators and the surviving time of integration, two of them (case a, b) are practically stable, so that the real system could be running in one of these configurations.

Via a new Hamiltonian expansion which is suitable for high-eccentricities planar three-body problem, we study the dynamics of the different configurations. We discuss the variations of eccentricities and resonant angles, explain the happenings of different evolving types, and give a criterion of the occurrence of the apsidal corotation.

The surfaces of section for the three types of motion are calculated and they reveal the stabilities of systems with or without the apsidal corotation. With these results we argue that the stability of the system is mainly due to the 3:1 MMR, and the apsidal corotation has only a limited contribution. We would like to mention that this method can also be applied to other extra-solar planetary systems with other mean motion resonances.

Numerical simulations suggest at least $`10`$ percent of systems are in a 3:1 MMR, while the Hamiltonian analyses give an upper limit of $`30\%`$, which however will drop down after considering the stability. We also list the initial conditions leading to this MMR.

Systems with different motion configurations have different energy ($`H`$) levels. The apsidal corotation happens when the Hamiltonian approaches the extreme value. So, if the two planets are captured into current configuration through orbital migration caused by the action of non-conservative forces, the system should have an extremum of energy so that the apsidal corotation happens \[Kley 2003\]. When the future observations would reveal more accurate properties of this system, we may consider what signatures of the migration are still presented in this system.

The masses of planets adopted in this paper are the values from the orbital solutions when assuming $`\mathrm{sin}i=1`$. As for the situations of $`\mathrm{sin}i<1`$, our initial analysis with the Hamiltonian get the very similar results in a wide range of $`\mathrm{sin}i`$. This is also consistent with the results in \[Beaugé et al 2003\]. We have also analysed the possible motion configurations and their stabilities if the initial eccentricities differ from the values adopted above. With the help of the Hamiltonian model, we find that a higher $`e_2`$ favors the establishment of a 3:1 MMR. Anyway, we’d like to leave more details of these to our future paper.

Last but not least, the general relativity effect may affects the secular dynamics of the 55 Cancri system, since the Companion B and C are quite close to the central star. For the innermost planet, the orbital precession caused by the general relativity effect, is calculated. Although this periastron shift, $`1.66\times 10^6`$ radians per orbit, is quite small, it’s about three times larger than that of Mercury in our Solar system.

## Acknowledgements

We thank Dr. Beaugé for helpful discussions. Appreciations also go to Dr. Tsevi Mazeh for the constructive comments and suggestions. This work is supported by Academy of Finland (Grant No. 44011). LYZ and YSS are also supported by the Natural Science Foundation of China (No. 10233020) and the Special Funds for Major State Basic Research Project (G200077303).


text:
# Winds and Outflows in Starburst Galaxies and AGN

## 1. Introduction

We present recent results on winds and outflows in starburst galaxies and active galactic nuclei (AGN) based on observations carried out with the Chandra X-ray observatory.

Outflows play a crucial role in enriching the nuclear environment with matter from the core region, carrying metals and possibly dust. They thus represent an important component when studying the recycling of interstellar and intergalactic matter.

In starburst galaxies, superwinds are driven by nuclear starbursts. Their study provides important information on the evolution of these winds, the metal content they carry, and the process of IGM enrichment with metals and/or dust.

In AGN, outflows may be driven by radiation pressure of the central continuum source. Components in outflow include possible accretion-disk winds and ionized absorbers. Chandra and XMM-Newton spectra of AGN with warm absorbers are rich in absorption features which allow to study the ionized material, its origin, evolution and interaction with the environment, in great detail (e.g., Simkin & Roychowdhury 2002).

## 2. Starburst-driven Superwinds, AGN environment: NGC 6240

### 2.1. The galaxy NGC 6240

Ultraluminous infrared galaxies (ULIRGs) are characterized by their huge luminosity output in the infrared, predominantly powered by super-starbursts and/or hidden AGN (e.g., Genzel et al. 1998, Sanders et al. 1988, 1999). Many distant SCUBA sources, massive and dusty galaxies, are believed to be ULIRG equivalents at high redshift. Local ULIRGs are therefore ideally suited to study the physics of galaxy formation and evolution (many ULIRGs are mergers), the processes of IGM enrichment, the physics of superwinds driven by the nuclear starbursts, and to asses the frequency of hidden AGN at their cores. The questions regarding the onset of starburst and AGN activity and their evolution in mergers are of fundamental importance for our understanding of AGN/black hole formation and evolution in general.

Here, we would like to concentrate on the galaxy NGC 6240, one of the nearest members of the class of (U)LIRGs, considered to be a key representative of its class. The galaxy is a merger in the process of forming an elliptical galaxy.

ROSAT observations (Schulz et al. 1998) showed that NGC 6240 is special in its highly luminous extended X-ray emission, making it one of the most luminous X-ray emitters in extended emission known among field galaxies (Komossa et al. 1998).

### 2.2. Chandra results

Given the complex nature of the X-ray emission of NGC 6240 with contributions from many components suggested from previous X-ray observations, spatially resolved spectroscopy is crucial to disentangle all contributing components, determine their nature, and derive their physical properties. Here, we report results from the first spatially resolved X-ray spectroscopy of NGC 6240, carried out with the ACIS-S instrument onboard the Chandra X-ray observatory (see also Komossa et al. 2003).

The Chandra image of NGC 6240 reveals a wealth of structure, changing in dependence of energy. A large part of the X-ray emission of NGC 6240 is extended, confirming previous results first obtained with the ROSAT instruments.

Basically, the extended emission of NGC 6240 shows three components: firstly, some weak, widely extended halo emission, which was already hinted for by our previous ROSAT observations, not discussed further here. Secondly, an emission component spatially coincident with the H$`\alpha `$ emission of NGC 6240. Thirdly, a more compact, harder, extended component, likely related to some nuclear starburst activity.

In the Chandra ACIS-S image, several X-ray ‘loops’ and knots are visible which correlate well with the H$`\alpha `$ emission of NGC 6240 (Fig. 2), and which are most prominent between 1 and 2.5 keV. Above 1.5 keV, X-ray emission from the direction of the northern nucleus of NGC 6240 starts to emerge. The hard X-ray image is dominated by emission from two compact sources, spatially coincident within the errors with the IR positions of the two nuclei of NGC 6240.

The extended X-ray emission of NGC 6240 of ”butterfly-like” shape and coincident with the H$`\alpha `$ emission, is generally well described by a MEKAL model with $`kT0.8`$ keV and absorption with a column density $`N_\mathrm{H}\mathrm{3\hspace{0.17em}10}^{21}`$ cm-2. Abundances are not well constrained, but fits typically yield $`Z0.1`$ solar when applying simple one-component thermal models.

In addition to the widely extended emission, a more compact nuclear component, still extended, is present. Representative of this emission, we have determined the spectrum of an emission-blob south-west of the Southern nucleus of NGC 6240. It is well fit by a MEKAL model of temperature $`kT3keV`$ and abundances of $`Z`$ solar.

The trend of decreasing abundances with increasing nuclear separation of the starburst component of Mrk 273 was noted by Xia et al. (2002). Whether this trend is true for NGC 6240, or whether the lower abundances inferred for the more extended starburst component of NGC 6240 are still an artifact of not applying more complex emission models (e.g., Komossa et al. 2000), has yet to be determined. It would require much deeper Chandra observation of NGC 6240 (and other (U)LIRGs) than presently available, in combination with more complex, possibly non-equilibrium emission models.

According to superwind models of Schulz et al. (1998), a mechanical input power of $`L_{\mathrm{mech}}=\mathrm{3\hspace{0.17em}10}^{43}`$ erg/s (derived from a SN-rate of 3/year) can, within 3 $`10^7`$ yrs, drive a single shell to an extent of $`R`$ 10 kpc within a medium of 0.1 cm-3 particle density. This model reproduces the observed X-ray luminosity of the extended, butter-fly shaped, starburst-related X-ray emission of NGC 6240 of $``$1042 erg/s.

The spectra of both X-ray nuclei of NGC 6240 are very hard in X-rays and show strong, $``$neutral iron K $`\alpha `$ lines. These properties identify both nuclei of NGC 6240 as active, and suggest a scattering geometry for both AGN in NGC 6240 (the core components are discussed in greater detail by Komossa et al. 2003). In addition to the 6.4 keV iron line, which is traced back to fluorescence from cold material in the AGN environment, a weaker line from highly ionized iron shows up, which could be related to supernova remnants or an ionized scatterer (see Komossa et al. 1998 for models of a warm scatterer in NGC 6240, and Boller et al. 2003 for recent XMM data which resolved the high-ionization iron-line complex in two components.)

## 3. Absorption and outflows in AGN

Cold and ‘warm’ absorption is ubiquitous in the AGN environment. X-ray absorption and emission features provide valuable diagnostics of the physical conditions in the X-ray gas and, in particular, allow to measure elemental abundances at high redshift, with profound consequences for our understanding of the star formation history in the early universe.

### 3.1. BAL quasars

Here, we would like to concentrate on BAL quasars. They are characterized by broad UV absorption lines. It has been suggested that these lines arise in a flow of gas which rises vertically from a narrow range of radii from the accretion disk. The flow then bends and forms a conical wind moving radially outwards (Elvis 2000). Variants of radiatively-driven disk-winds were explored (e.g., Proga et al. 2000, and references therein). In some of these models, an X-ray absorber shields the wind downstream from soft X-rays, allowing resonant-line driving to remain effective and accelerate the outflowing BAL wind up to $``$0.1c.

These high-velocity BAL winds may play a crucial role in metal-enrichment of the interstellar medium of the galaxy, possibly even the IGM.

### 3.2. The quasar APM 08279+5255

The high-redshift BAL quasar APM 08279+5255 is magnified by a gravitational lense and is among the most luminous objects in the universe (even after correction for lensing). It has one of the best-measured X-ray spectra among BALs. The recent XMM-Newton observation led to the detection of a strong absorption feature of ionized iron, interpreted as K-edge, arising from a warm absorber of high column density (Hasinger et al. 2002), and variable on short timescales, when compared with the Chandra observation (Chartas et al. 2002).

We concentrate here on the aspect of deriving metal abundances at high redshift (for an application to measure cosmological parameters, see e.g. Alcaniz et al. 2003). Using our best-fit warm absorber model of the XMM spectrum of APM 08279+5255, we infer an abundance of Fe/O $``$ 3 $`\times `$ solar which is rather high, given the high redshift of the object.

The element iron plays a special role in chemical enrichment scenarios, because its production is delayed relative to other elements (e.g., Fig 1 of Hamann & Ferland 1993), since it is believed to be mostly produced in supernovae of type Ia. Taken at face value, the XMM observation of APM 08279+5255 imply an efficient iron production mechanism in the early universe. The unusual strength of the iron features and the indications for strong variability make APM 08279+5255 an excellent target for simultaneous, deep follow-up observations with Chandra and XMM-Newton.

While Fe X-ray features have only been observed in very few BALs so far, with future X-ray missions like XEUS (Arnaud et al. 1998) we will be able to observe more objects with high accuracy. The iron absorption lines and K-edges provide a unique probe of matter at high redshift because, firstly, they are easy to measure even, or, particularly, at high redshift $`z`$, and secondly, Fe(/O) abundance measurements in the early universe are important for constraints on the early star formation history.

Science issues that we will be able to address in detail with future X-ray missions for the first time, particularly at high redshift, include the determination of metal abundances of X-ray absorbers by detection of metal absorption edges, the determination of the composition of dust mixed with cold and ionized gas (K-edges of metals in cold dust and cold gas will be resolvable from each other for the first time), measurements of the velocity field of the gas, and the utilization of these results to investigate the effects of reprocessing of gas and dust in AGN from high to low redshift.

#### Acknowledgments.

I, St.K., am grateful to Hartmut Schulz for intense collaborations on NGC 6240 ever since I started in astrophysics, for sharing his wide knowledge, and for ongoing discussions on many topics of astrophysics. Hartmut passed away in August 2003, after a short, severe illness.
We thank R. van der Marel for providing the HST H$`\alpha `$ image of NGC 6240.

## References

Alcaniz J.S., Lima, J.A.S., Cunha, J.V. 2003, MNRAS, 340, L39

Arnaud, M., et al. 1999, The XEUS Science Case, ESA

Boller, T., et al. 2003, astro-ph/0307326

Chartas, G., Brandt, W.N., Gallagher, S.C., Garmire, G.P., 2002, ApJ, 597, 169

Elvis, M., 2000, ApJ, 545, 63

Genzel, R., Lutz, D., Sturm, E., et al. 1998, ApJ, 498, 579

Gerssen, J., van der Marel, A.P., Axon, D., Mihos, C., Hernquist, L., Barnes, J.E., 2001, ASP Conf. Ser., 249, 665

Hamann, F., Ferland, G.J. 1993, ApJ, 418, 11

Hasinger, G., Schartel, N., Komossa, S. 2002, ApJ, 573, L80

Komossa, S., Schulz, H., Greiner, J. 1998, A&A, 334, 110

Komossa, S., Breitschwerdt, D., Böhringer, H., Meerschweinchen, J. 2000, Ap&SS, 272, 295

Komossa, S., Burwitz, V., Hasinger, G., Predehl, P., Kaastra, J.S., Ikebe, Y., 2003, ApJ, 582, L15

Proga, D., Stone, J.M., Kallman, T.R., 2000, ApJ, 543, 686

Sanders, D.B., Soifer, B.T., Elias, J.H., et al. 1988, ApJ, 325, 74

Sanders, D.B. 1999, Ap&SS, 266, 331

Schulz, H., Komossa, S., Berghöfer T., Boer, B. 1998, A&A, 330, 823

Simkin, M.V., Roychowdhury, V.P., 2002, cond-mat/212043

Xia, X.Y., Xue, S.J., Mao, S., et al. 2002, ApJ, 564, 196