text:
# Representative Pathways of Excitation Migration in Photosystem I

## Introduction

Photosynthesis is carried out by pigment-protein complexes embedded in cell membranes. An aggregate of interacting pigments held in a fixed arrangement in such a complex absorbs light and delivers the resulting electronic excitation to a reaction center, which uses the excitation energy to produce a charge separation across the cell membrane Blankenship (2002). The transmembrane potential induced by this charge separation is later used for the synthesis of ATP.

Among the photosynthetic complexes, the photosynthetic unit of purple bacteria has been most extensively studied; for a review see Hu et al. (2002). Recently, a high-resolution structure of photosystem I (PSI) has been obtained from cyanobacterium Synechococcus elongatus Jordan et al. (2001). PSI along with photosystem II constitutes the machinery of oxygenic photosynthesis in plants, green algae, and cyanobacteria. PSI contains an aggregate of 96 chlorophylls, including the special pair P700, where charge separation is initiated. The chlorophylls are arranged without an apparent order, except for an apparent pseudo-$`C_2`$ symmetry, which is in contrast to the highly symmetric arrangement of bacteriochlorophylls in the photosynthetic unit of purple bacteria. The rather random arrangement of chlorophylls in PSI poses the question which pathways excitation migration follows towards P700 after any of its chlorophylls has absorbed light.

The rates of inter-chlorophyll excitation transfer in PSI have been calculated based on Förster theory Byrdin et al. (2002); Şener et al. (2002); Damjanović et al. (submitted). Excitation follows a stochastic trajectory along the excitation transfer network given by these rates. However, the obscure pattern of the excitation transfer network in PSI (cf. Fig. 6 in Ref. Byrdin et al. (2002)) instills a need for a simpler and more distilled picture of excitation migration. Also, the aspect that the special pair P700 is the target of the excitation migration process is not incorporated in the excitation transfer network solely given by the transfer rates. Therefore, we seek paths directed towards P700 that are representative of all events of excitation migration towards P700 in a way that actual migration trajectories can be considered as noisy trajectories around those representative paths.

The task of finding such representative pathways is also encountered in the so-called reaction-path problem Elber (1996); Straub (2001). Recently, a new approach to the reaction-path problem was suggested Park and Schulten (submitted), where mean first-passage times were employed as reaction coordinates and reaction paths were constructed from them. The suggested method is applicable to discrete systems such as excitation migration in photosynthetic complexes. In this letter we apply the method to find excitation migration pathways in PSI.

## Theory

Excitation migration in PSI is a stochastic process governed by the rates of inter-chlorophyll transfer, dissipation (internal conversion of excitation to heat), and the charge separation at P700. We build a master equation describing the process, calculate the mean first-passage time for excitation migration from each chlorophyll to P700, and then find pathways representative of the excitation migration in the excitation transfer network of the chlorophyll aggregate towards the P700 pair.

The process of excitation migration can be described in terms of the probability $`P_{ij}(t)`$ that an excitation, initiated by light absorption at chlorophyll $`j`$, resides at chlorophyll $`i`$ after time interval $`t`$. The inter-chlorophyll transfer rates $`T_{ij}`$ from chlorophyll $`i`$ to chlorophyll $`j`$ are calculated as explained in Ref. Şener et al. (2002). For the calculation of $`T_{ij}`$, we use the recently obtained site energies for chlorophylls Damjanović et al. (submitted) and the inter-chlorophyll electronic couplings given by a full Coulomb computation that includes all multipole contributions to the coupling Şener et al. (2002). The dissipation rate is assumed to be the same ($`k_{\mathrm{diss}}=(1\mathrm{ns})^1`$) at all the chlorophylls. Since we are here interested in paths along which excitation migrates towards P700, it is sufficient to consider only first passages of excitation to P700. (The first passage corresponds to the first term in the so-called sojourn expansion whose higher-order terms cover subsequent escapes and returns to P700 Şener et al. (2002).) Thus, the charge-separation rate is not needed in the model, and the transfer rates from P700 to the other chlorophylls can be set to zero. Collecting these rates, we obtain the master equation Ritz et al. (2001)

$$\frac{d}{dt}P_{ij}(t)=_kK_{ik}P_{kj}(t),$$
(1)
where the $`96\times 96`$ transition matrix $`K_{ij}`$ is

$$K_{ij}=\{\begin{array}{cc}\hfill T_{ji}\delta _{ij}_kT_{jk}\delta _{ij}k_{\mathrm{diss}}\mathrm{for}& j\mathrm{P700}\hfill \\ \hfill 0\mathrm{for}& j\mathrm{P700}\hfill \end{array}.$$
(2)
From the master equation (Eq. 1) we derive a backward master equation which is used for calculating the mean first-passage time Nadler and Schulten (1985); Gardiner (1985). Since excitation migration is a Markov process, the probability $`P_{ij}(t)`$ can be expressed as

$$P_{ij}(t)=_kP_{ik}(s)P_{kj}(ts)$$
(3)
for an arbitrary intermediate time $`s`$ ($`0<s<t`$). Taking the derivative with respect to $`t`$ leads to

$`{\displaystyle \frac{d}{dt}}P_{ij}(t)`$ $`=`$ $`{\displaystyle _k}P_{ik}(s){\displaystyle \frac{d}{dt}}P_{kj}(ts)`$ (4)
$`=`$ $`{\displaystyle _k}P_{ik}(s){\displaystyle _l}K_{kl}P_{lj}(ts).`$
We now take the limit of $`s(t0)`$ and obtain the backward master equation

$$\frac{d}{dt}P_{ij}(t)=_kP_{ik}(t)K_{kj},$$
(5)
where we have used $`lim_{s(t0)}P_{ij}(ts)=\delta _{ij}`$.

In the model constructed above, excitation can exit the system by two means: the delivery to P700 and dissipation. We emphasize again that we only consider first passages to P700. Thus, it is convenient to consider a subsystem of 94 chlorophylls, excluding the P700 pair of chlorophylls. The subsystem is described by the master equation

$$\frac{d}{dt}P_{\alpha \beta }(t)=_\gamma K_{\alpha \gamma }P_{\gamma \beta }(t)$$
(6)
or the backward master equation

$$\frac{d}{dt}P_{\alpha \beta }(t)=_\gamma P_{\alpha \gamma }(t)K_{\gamma \beta },$$
(7)
where Greek subscripts indicate that the P700 chlorophylls are not included. The $`94\times 94`$ transition matrix $`K_{\alpha \beta }`$ is obtained from the original $`96\times 96`$ transition matrix $`K_{ij}`$ by eliminating the rows and columns belonging to the P700 chlorophylls.

Consider an excitation initiated at chlorophyll $`\alpha `$ at time zero. The probability that the excitation reaches P700 between time $`t`$ and $`t+dt`$ is equal to $`_\beta dt\xi _\beta P_{\beta \alpha }(t)`$, where $`\xi _\beta =_{i\mathrm{P700}}T_{\beta i}`$ is the transfer rate to P700 from chlorophyll $`\beta `$. Thus, the (conditional) mean first-passage time $`\tau _\alpha `$, namely the average exit time given that the exit has occurred through a delivery to P700 (not through dissipation), is given as

$$\tau _\alpha =_0^{\mathrm{}}𝑑tt_\beta \xi _\beta P_{\beta \alpha }(t)/\varphi _\alpha $$
(8)
$$\varphi _\alpha =_0^{\mathrm{}}𝑑t_\beta \xi _\beta P_{\beta \alpha }(t).$$
(9)
Here $`\varphi _\alpha `$ is the total probability that the excitation is eventually delivered to P700. Since the excitation may be dissipated before it reaches P700, $`\varphi _\alpha `$ is less than one and the denominator in Eq. 8 is necessary for normalization. From Eqs. 78, and 9, we obtain equations for $`\varphi _\alpha `$ and $`\tau _\alpha `$:

$$_\alpha \varphi _\alpha K_{\alpha \beta }=\xi _\beta $$
(10)
$$_\alpha \tau _\alpha \varphi _\alpha K_{\alpha \beta }=\varphi _\beta .$$
(11)
The solutions to these equations are given in terms of the inverse matrix $`K_{\alpha \beta }^1`$ of the matrix $`K_{\alpha \beta }`$:

$$\varphi _\alpha =_\beta \xi _\beta K_{\beta \alpha }^1$$
(12)
$$\tau _\alpha =_\beta \varphi _\beta K_{\beta \alpha }^1/\varphi _\alpha .$$
(13)
The mean first-passage time $`\tau _i`$ ($`\tau _i=0`$ for $`i\mathrm{P700}`$) is a natural measure of the progress in the excitation migration towards the P700 pair Park and Schulten (submitted); chlorophylls associated with shorter mean first-passage times are closer to the target (P700) in the migration process. As suggested in Ref. Park and Schulten (submitted), we construct excitation migration paths according to the scheme that a path going through chlorophyll $`i`$ chooses the next chlorophyll $`j`$ such that $`T_{ij}(\tau _i\tau _j)`$ is maximized. The resulting paths can be considered representative because the mean first-passage time $`\tau _i`$ is an average over all events of the excitation migration from chlorophyll $`i`$ to P700. For a transfer step from chlorophyll $`i`$ to chlorophyll $`j`$, the quantity $`1/T_{ij}`$ may be interpreted as a cost, and $`\tau _i\tau _j`$ as a gain. The scheme then amounts to maximizing the gain-cost ratio for each step.

## Results and Discussion

The above method yields the mean first-passage time and the paths shown in Fig. 1. The 96 chlorophylls are sorted in the order of increasing mean first-passage time and are listed in Fig. 1(c). The excitation migration paths shown in Fig. 1(a) exhibits a network without an apparent order, as expected from the arrangement of the chlorophylls. There are various paths, including six direct paths to P700 (from chlorophylls eC-A2, eC-B2, A24, A26, B5, and B24) and the most complicated path composed of nine steps (A21$``$PL1$``$A22$``$A30$``$B37$``$B38$``$B39$``$eC-B3$``$eC-A2$``$eC-B1). We here emphasize again that excitation does not necessarily follow these paths, but that they should be understood as representative paths.

As can be seen from Fig. 1(b), most of the chlorophylls (9th to 96th) have mean first-passage times around or above 10 ps, with the average over all being 12.5 ps. Chlorophylls A40 and B39 lie between these peripheral chlorophylls and the reaction-center chlorophylls which coincide with the six chlorophylls with the shortest mean first-passage times. This supports to a certain extent the suggestion that chlorophylls A40 and B39 connect the reaction-center chlorophylls and the peripheral chlorophylls Jordan et al. (2001). But, among the found excitation migration paths to P700, many go through neither chlorophyll A40 nor B39. Therefore, these two chlorophylls are not the only connection between the reaction center and the periphery, and should not be considered as bottlenecks.

The last chlorophyll in Fig. 1(c), chlorophyll M1, raises a question as its mean first-passage time (18.8 ps) is much longer than the second longest one (15.6 ps). This does not imply that excitation is trapped at this chlorophyll; 18.8 ps is still much shorter than the dissipation time of 1 ns. But, it seems inappropriate that one chlorophyll is located relatively farther from the rest. PSI exists as a trimer in vivo Boekema et al. (1987); Jordan et al. (2001) and chlorophyll M1 lies close to the boundary between monomers. By computing the electronic couplings between chlorophyll M1 and the chlorophylls in the next monomer, we find that chlorophyll M1 is coupled to chlorophyll A30 in the next monomer with the coupling of 52.9 cm-1, which is much stronger than the strongest coupling it has within its own monomer (6.7 cm-1 with B8). Chlorophyll M1 functionally belongs to the next monomer, not its own monomer.

In PSI, it is experimentally known that seven to eleven chlorophylls absorb light at longer wavelengths than the special pair P700 Pålsson et al. (1998); Zazubovich et al. (2002). These ‘red chlorophylls’ have not been identified with certainty, but some candidates (marked with red dots in Figs. 1(a) and (c)) have been proposed Jordan et al. (2001); Byrdin et al. (2002); Damjanović et al. (submitted); Şener et al. (2002). Most of the candidates are found to be located rather close to P700 in terms of not only spatial distance but also the mean first-passage time. Chlorophylls A38, A39, B37, and B38 have the shortest mean first-passage times among peripheral chlorophylls (except for B39 and A40). With regard to the representative paths shown in Fig. 1(a), chlorophylls A24, A32, A38, B22, B24, and B26 are one step away from the reaction-center chlorophylls, and chlorophylls A39, B6, and B38 are one step away from the ‘connecting chlorophylls’, A40 and B39.

In conclusion, we have found representative paths of excitation migration in PSI based on mean first-passage times. The mean first-passage time and the paths provide a complete yet distilled picture of the excitation migration towards P700 as illustrated in Fig. 1. We expect that our methodology will be useful for various light-harvesting complexes as more of their high-resolution structures become available.

## Acknowledgments

This work has been supported by National Institute of Health grant PHS 5 P41 RR05969.


format: text

text:
# Comment on “Supersymmetry in the half-oscillator - revisited”

## Abstract

We point out the flaw in the analysis of Gangopadhyaya and Mallow, hep-th/0206133, where it is claimed that supersymmetry is broken in the SUSY half-oscillator, even with a regularization respecting supersymmetry.

In an earlier paper das , we had shown that supersymmetry, in quantum mechanical theories with singular potentials or nontrivial boundaries, is preserved if the system is regularized in a manner respecting supersymmetry. This is in contrast to the earlier claims jevicki ; roy ; cooper  that supersymmetry is broken in such systems. We had shown, in particular, that if the superpotential is regularized, as opposed to the conventional wisdom of regularizing the potential, supersymmetry is maintained in such systems in a natural manner. The reason for this is quite clear. A regularized superpotential leads to a pair of supersymmetric Hamiltonians for every value of the regularization parameter whereas a conventionally regularized potential does not lead to a pair of supersymmetric Hamiltonians for any value of the regularization parameter.

In a recent paper, Gangopadhyaya and Mallow asim  claim that of the examples treated in das , the supersymmetric half-oscillator does not possess supersymmetry in spite of using the supersymmetric regularizations proposed in das . In fact, their analysis is exactly the same as given in das . However, their final solution of the boundary condition and, therefore, their conclusion is incorrect for a very simple reason that we wish to point out.

Let us recall that the SUSY half-oscillator was first studied in roy  where the authors concluded that supersymmetry is broken. Physically, their result can be understood as follows. If we require the wavefunctions to vanish at the origin (as these authors do and as would be the case when one conventionally regularizes the potential), only odd solutions of the harmonic oscillator are allowed. On the other hand, supersymmetry requires the partner wavefunctions to have opposite parity and since this is not allowed by the boundary conditions (which follows from the particular regularization used), supersymmetry must be broken. The problem with this conclusion, as pointed out in das , is that such a treatment of the system corresponds to regularizing the potential in a manner which itself breaks supersymmetry and, therefore, the result of such an analysis can be an artifact of the regularization used. In stead, we had proposed to regularize the superpotential directly as

$$W(x)=\omega x\theta (x)+c\theta (x)$$
(1)
where $`c`$ is the regularization parameter and we take $`c+\mathrm{}`$ at the end. This superpotential leads to the regularized supersymmetric potentials:

$`V_+(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[(\omega ^2x^2\omega )\theta (x)+c^2\theta (x)c\delta (x)]`$
$`V_{}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[(\omega ^2x^2+\omega )\theta (x)+c^2\theta (x)+c\delta (x)]`$ (2)
The crucial difference between this and the earlier treatment of this problem roy  is the appearance of Dirac delta functions in the potentials. The delta function is attractive for one of the Hamiltonians while it is repulsive for the other. It is clear, therefore, that, in the presence of an attractive delta function potential, the usual argument for the non existence of even solutions is not automatic. In fact, we had shown explicitly that supersymmetry is maintained in the limit $`c+\mathrm{}`$ in this theory. We had also proposed an alternate regularization in das  that avoids the appearance of delta functions, and had shown that the end result is the same if the regularization maintains supersymmetry.

It is surprising, therefore, that Gangopadhyaya and Mallow asim , who carry out exactly our analysis would arrive at a different conclusion, namely, that supersymmetry is broken in this system. Let us note that the boundary condition (eq. (23) in das  and eq. (10) in asim ) is given by (we choose $`\omega =\frac{1}{2}`$ for simplicity)

$$\sqrt{2}\frac{\mathrm{\Gamma }(\frac{1}{2}ϵ)}{\mathrm{\Gamma }(ϵ)}=\sqrt{c^22ϵ}c$$
(3)
where $`ϵ`$ represents the energy of the eigenstates. Gangopadhyaya and Mallow claim that $`ϵ=0`$ is not a solution of the above equation, in which case supersymmetry will be broken. It is, of course, obvious that $`ϵ=0`$ is a solution of (3) for any value of $`c`$ and, in particular, for $`c+\mathrm{}`$, but let us look at this more carefully to understand the flaw in the analysis of asim .

Let us define, for simplicity,

$$F(ϵ)\sqrt{2}\frac{\mathrm{\Gamma }(\frac{1}{2}ϵ)}{\mathrm{\Gamma }(ϵ)}\mathrm{and}G(ϵ,c)\sqrt{c^22ϵ}c,$$
(4)
so that we can represent the boundary condition (3) as

$$F(ϵ)=G(ϵ,c)$$
(5)
In a quantum theory, the boundary condition (5) may not hold for all values of $`ϵ`$. The fact that energy is quantized means that it can be satisfied only for some discrete values of $`ϵ=ϵ_n`$. Let us, in particular, investigate whether $`ϵ=0`$ is a solution of the equation, and if it is, whether it remains a solution in the limit $`c+\mathrm{}`$.

Let us note that both $`F(ϵ)`$ and $`G(ϵ,c)`$, as complex functions of $`ϵ`$, are analytic around $`ϵ=0`$. For $`F(ϵ)`$, we can see this from the fact that the function $`\frac{1}{\mathrm{\Gamma }(ϵ)}`$ is entire in the whole $`ϵ`$ plane with simple zeroes at the points $`ϵ=n`$, $`(n=0,1,2,\mathrm{})`$, and the function $`\mathrm{\Gamma }(\frac{1}{2}ϵ)`$ is also analytic around $`ϵ=0`$, with $`\mathrm{\Gamma }(\frac{1}{2})=\sqrt{\pi }`$. For $`G(ϵ,c)=\sqrt{c^22ϵ}c`$, the same is also true.

Since both these functions are analytic around $`ϵ=0`$, a Taylor expansion around this value leads to,

$`F(ϵ)`$ $``$ $`0+\sqrt{2\pi }ϵ+𝒪(ϵ^2)`$ (6)
$`G(e,c)`$ $``$ $`0{\displaystyle \frac{1}{c}}ϵ+𝒪(ϵ^2)`$ (7)
The fact that $`ϵ=0`$ is a solution of equation (5) is a consequence of the fact that the zeroth order terms in the Taylor expansion in equations (6) and (7) coincide for that value. Note that this result is independent of $`c`$ and, therefore, also holds in the limit $`c+\mathrm{}`$.

The mistake in the analysis of Gangopadhyaya and Mallow lies in the fact that, for whatever reason, they impose that the coefficient of the linear terms in the Taylor expansion should match, in addition to the zeroth order terms. Indeed, the equation below equation (11) (for $`\omega =\frac{1}{2}`$) in asim  is precisely the requirement of the equality for the linear terms in equations (6) and (7). This imposition, however, is unjustified, since that would tantamount to requiring

$$F(0)=G(0,c),\mathrm{and}F^{}(0)=G^{}(0,c)$$
(8)
The boundary condition (5) does not have to hold for all values of $`ϵ`$, not even in a neighborhood of $`ϵ=0`$. It should hold only for discrete values of $`ϵ`$. For the case $`ϵ=0`$, this translates into the fact that only the zeroth order terms in the Taylor expansion of $`F`$ and $`G`$ should be equal, but no extra conditions must be imposed on higher order terms. The incorrect conclusion in asim  results from the imposition of this additional condition. Supersymmetry is, in fact, unbroken in this theory, when analyzed properly.

Acknowledgment:

This work was supported in part by US DOE Grant number DE-FG 02-91ER40685.


text:
# Polarized Parton Distributions in the Valon Model

## 1 INTRODUCTION

The idea of quark cluster is not new. In valon model, the hadron is envisaged as a bound state of valence quark cluster called valon . For example, the bound state of proton consists of two “up” and one “down” valon. These valons thus bear the quantum numbers of respective valence quarks. Hwa  found evidence in the deep inelastic neutrino scattering data that suggested their existence and applied it to a variety of phenomena. Hwa  had also successfully formulated a treatment of the low-$`p_T`$reactions based on the structural analysis of the valons. Arash  has recently applied the valon model to study the contributions of the constituents of the proton to its spin in deep inelastic scattering (DIS). Results of the model calculation agree well with the European Muon Collaboration (EMC) results . Arash adopt the valon point of view and calculate various contributions to the spin of the proton.

The valon description of the nucleon was developed by Hwa  some years ago and was successfully used  to describe the contribution of constituents of the proton to its spin in order to understand the so-called spin crisis. In this paper, we investigate the calculation of polarized parton density of the nucleon using an improved model and apply the model to the study of the spin of nucleon. In Sec. 2, polarized valon distributions are derived. In this section we calculate polarized valon distributions as a function of $`y`$ from unpolarized valon distributions. The moment of polarized parton distributions and using inverse Mellin transformation are discussed in Sec.3. In Sec. 4, we examine the polarized structure function, $`xg_1`$, for proton, neutron and deuteron and compare with the experimental data. Also our prediction for the first moment of these structure functions presents in this section too.

## 2 UNPOLARIZED AND POLARIZED VALON DISTRIBUTIONS

The polarized valon is defined to be a dressed polarized valence quark in QCD with the cloud of polarized gluons and sea quarks which can be resolved by high $`Q^2`$ probes. In a scattering process the virtual emission and absorption of gluons in a valon become bremsstrahlung and pair creation. The polarized structure function of a valon is determined by gluon bremsstrahlung and pair creation in QCD. At sufficiently low value of $`Q^2`$ the internal polarized structure of a valon cannot be resolved and hence it behaves as a polarized valence quark.

To facilitate the phenomenological analysis Hwa  assumed a simple form for the exclusive valon distribution in unpolarized proton which is

$$G_{UUD/p}(y_1,y_2,y_3)=g_p(y_1y_2)^\alpha y_3^\beta \delta (y_1+y_2+y_31)$$
(1)
$`y_i`$ is the momentum fraction of the $`i`$th valon , and $`g_p`$ is

$$g_p=[B(\alpha +1,\beta +1)B(\alpha +1,\alpha +\beta +2)]^1$$

where $`B(m,n)`$ is the beta function. The valon distributions are obtained by integration

$`G_{U/p}(y)`$ $`=`$ $`{\displaystyle 𝑑y_2𝑑y_3G_{UUD/p}(y,y_2,y_3)}`$
$`=`$ $`g_pB(\alpha +1,\beta +1)y^\alpha (1y)^{\alpha +\beta +1}`$
$`G_{D/p}(y)`$ $`=`$ $`{\displaystyle 𝑑y_1𝑑y_2G_{UUD/p}(y_1,y_2,y)}`$
$`=`$ $`g_pB(\alpha +1,\alpha +1)y^\beta (1y)^{2\alpha +1}`$
Recently, R.C. Hwa  have recalculated the unpolarized valon distribution in proton with a new set of parameters. The new values of $`\alpha `$ , $`\beta `$ are found to be $`\alpha =1.76`$ and $`\beta =1.05`$. An alternative analysis of polarized structure function has been performed by Gehrmann and Stirling. This was done in the same spirit as the unpolarized analysis , namely the polarized valon distribution chose to be of the general form

$$\delta G_{U/p}(y)=\delta F_U(y)\times G_{U/p}(y)$$

$$\delta G_{D/p}(y)=\delta F_D(y)\times G_{D/p}(y)$$
(4)
where $`\delta ϝ_j(y)`$ is define:

$$\delta F_j(y)=N_jy^{\alpha _j}(1y)^{\beta _j}(1+\gamma _jy+\delta _jy^{0.5})$$
(5)
where the subscript $`j`$ refer to $`U,D`$, and $`G_{\frac{v}{p}}(y)`$ is defined in Eqs. (2,3). In order to write down the polarized valon distribution in a proton we assume that at some $`Q=Q_0`$ a valon behaves as a polarized valence quark. The polarized quark distribution in a proton at a starting scale of $`Q_0^2=0.26GeV^2`$ is parameterized by GRV group . We take that parameterization and extract all of the parameters as defined in Eq.(5). The corresponding parameters are given in Table I. In Figure (1) we present unpolarized and polarized valon distributions as a function of $`y`$.

## 3 POLARIZED PARTON DISTRIBUTIONS

We now go to the $`n`$-moment space and define the moment of polarized valon distribution $`\mathrm{\Delta }M_{v/p}(n)`$ as:

$$\mathrm{\Delta }M_{v/p}(n)_0^1y^{n1}\delta G_{v/p}(y)𝑑y$$
(6)
the subscription $`v`$ stands for kind of $`U,D`$ polarized valon. For the moments of polarized singlet and non singlet distributions we shall use, the leading order solutions of the renormalization group equation in QCD. They can be expressed entirely in terms of the evolution parameter $`s`$

$$s=\mathrm{ln}\frac{\mathrm{ln}Q^2/\mathrm{\Lambda }^2}{\mathrm{ln}Q_0^2/\mathrm{\Lambda }^2}$$
(7)
where $`Q_0,\mathrm{\Lambda }`$ are scale parameters to be determined by experiments. From the theoretical standpoint both $`Q_0`$ and $`\mathrm{\Lambda }`$, depend on the order of the moments $`n.`$ In this work we have assumed that they are independent of n and we fixed it by $`Q_0^2=0.26GeV^2`$ and $`\mathrm{\Lambda }=0.232GeV`$.

The moments of singlet and non singlet polarized are

$$\mathrm{\Delta }M^{NS}(n,Q^2)=\mathrm{exp}(\delta d_{NS}^{(0)n}s)$$
(8)
$$\mathrm{\Delta }M^S(n,Q^2)=\frac{1}{2}(1+\delta \rho )\mathrm{exp}(\delta d_+s)+\frac{1}{2}(1\delta \rho )\mathrm{exp}(\delta d_{}s)$$
(9)
where $`\delta \rho `$ and other associated parameters are as follows:

$`\delta \rho `$ $`=`$ $`{\displaystyle \frac{\delta d_{NS}^{(0)n}\delta d_{gg}^{(0)n}+4f\delta d_{qg}^{(0)n}}{\mathrm{\Delta }}}`$
$`\mathrm{\Delta }`$ $`=`$ $`\delta d_+\delta d_{}=\sqrt{(\delta d_{NS}^{(0)n}\delta d_{gg}^{(0)n})^2+8f\delta d_{qg}^{(0)n}\delta d_{gq}^{(0)n}}`$
$`\delta d_\pm `$ $`=`$ $`{\displaystyle \frac{1}{2}}(\delta d_{NS}^{(0)n}+\delta d_{gg}^{(0)n}\pm \mathrm{\Delta })`$ (10)
$`b`$ $`=`$ $`{\displaystyle \frac{332f}{12\pi }}`$
and also the anomalous dimensions $`\delta d_{ij}^{(0)n}`$ are simply the n-th moment of polarized LO splitting function are given by  and $`f`$ is the number of active flavors.

$$\delta d_{NS}^{(0)n}=\frac{1}{3\pi b}[3\frac{2}{n(n+1)}+4S_1(n)]$$

$$\delta d_{qg}^{(0)n}=\frac{1}{4\pi b}\frac{n1}{n(n+1)}$$

$$\delta d_{gq}^{(0)n}=\frac{2}{\pi b}\frac{1}{3}\frac{n+2}{n(n+1)}$$

$$\delta d_{gg}^{(0)n}=\frac{1}{2\pi b}\{3[\frac{11}{6}+\frac{4}{n(n+1)}2S_1(n)]\frac{2}{3}\frac{f}{2}\}.$$
(11)
Here $`S_1(n)=_{j=1}^n\frac{1}{j}=\psi (n+1)+\gamma _E,\psi (n)=\frac{\mathrm{\Gamma }^{}(n)}{\mathrm{\Gamma }(n)}`$ and $`\gamma _E=0.577216.`$

The moment of polarized $`u`$ and $`d`$-valence quark in a proton are:

$$\mathrm{\Delta }M_{u_v}(n,s)=2\mathrm{\Delta }M^{NS}(n,s)\mathrm{\Delta }M_{U/p}$$

$$\mathrm{\Delta }M_{d_v}(n,s)=\mathrm{\Delta }M^{NS}(n,s)\mathrm{\Delta }M_{D/p}$$
(12)
and the moment of polarized $`\mathrm{\Sigma }`$ and gluon distributions are as follows:

$$\mathrm{\Delta }M_\mathrm{\Sigma }(n,s)=\mathrm{\Delta }M^s(n,s)(2\mathrm{\Delta }M_{U/p}+\mathrm{\Delta }M_{D/p})$$

$$\mathrm{\Delta }M_g(n,s)=\mathrm{\Delta }M_{gq}(n,s)(2\mathrm{\Delta }M_{U/p}+\mathrm{\Delta }M_{D/p})$$
(13)
where $`\mathrm{\Delta }M_{gq}(n,s)`$ is the quark-to-gluon evolution function give by:

$`\mathrm{\Delta }M_{gq}(n,s)`$ $`=`$ $`{\displaystyle \frac{(\delta d_+\delta d_{NS}^{(0)n})}{\mathrm{\Delta }}}(1{\displaystyle \frac{\delta d_{}\delta d_{NS}^{(0)n}}{2f\delta d_{qg}^{(0)n}}})\mathrm{exp}d_+s`$ (14)
$`{\displaystyle \frac{(\delta d_{}\delta d_{NS}^{(0)n})}{\mathrm{\Delta }}}(1{\displaystyle \frac{\delta d_+\delta d_{NS}^{(0)n}}{2f\delta d_{qg}^{(0)n}}})\mathrm{exp}d_{}s`$
In Figure (2) we present these moments as a function of $`n`$ for some values of $`Q^2`$. To obtain the $`x`$dependence of structure functions and parton distributions, usually required for practical purposes, from the above $`n`$dependent exact analytical solutions in Mellin-moment space, one has to perform a numerical integral in order to invert the Mellin-transformation in (8) according to

$$\delta f(x,Q^2)=\frac{1}{\pi }_0^{\mathrm{}}𝑑zIm[e^{i\varphi }x^{cze^{i\varphi }}\mathrm{\Delta }f^{n=c+ze^{i\varphi }}(Q^2)]$$
(15)
where the contour of integration, and thus the value of $`c`$, has to lie to the right of all singularities of $`\mathrm{\Delta }f^n(Q^2)`$ in the complex $`n`$-plane, i.e., $`c0`$ since according to Eq. (13) the dominant pole of all $`\delta d_{ij}^n`$ is located at $`n=0`$. For all practical purpose one may choose $`c1,\varphi =135^{}`$ and an upper limit of integration, for any $`Q^2`$ , of about $`5+10/\mathrm{ln}x^1`$, instead of $`\mathrm{}`$, which suffices guarantee and stable numerical results . In figure (3) our predictions for $`x`$dependence spin densities for quarks and gluons at $`Q^2=2,3,5,10`$ $`GeV^2`$ in a proton are shown.

Furthermore, in the polarized parton distribution we assume an unbroken $`SU(3)_f`$ symmetric sea,

$$\delta \overline{q}(x,Q^2)\delta \overline{u}=\delta u_{sea}=\delta \overline{d}=\delta d_{sea}=\delta \overline{s}=\delta s_{sea}$$
(16)
Note that the proton include partons with positive and negative helicity which carry of fraction $`x`$ of the proton momentum and to whom one can associate quark densities $`q_+(x,Q^2)`$ and $`q_{}(x,Q^2)`$. The difference

$$\delta q(x,Q^2)=q_+(x,Q^2)q_{}(x,Q^2)$$
(17)
measures how much the parton of flavor $`q`$ “remembers” of its parent parton polarization. Similarly, one may define

$$\delta \overline{q}(x,Q^2)=\overline{q}_+(x,Q^2)\overline{q}_{}(x,Q^2)$$
(18)
for antiquarks.

## 4 POLARIZED STRUCTURE FUNCTION, $`𝐠_1`$

We used polarized parton distribution from valon model for extracting structure function $`g_1`$. Within this model according Sec. (3), we calculated all of the polarized parton distributions in a proton. Now it is very easy to used these results for calculating $`g_1^p`$. In the quark model, to leading order (LO) in QCD, $`g_1^p`$ can be written as a linear combination of $`\delta q`$ and $`\delta \overline{q}`$ ,

$$g_1^p(x,Q^2)=\frac{1}{2}\underset{q}{}e_q^2[\delta q(x,Q^2)+\delta \overline{q}(x,Q^2)]$$
(19)
where $`e_q`$ are the electric charges of the (light) quark-flavors $`q=u,d,s`$.

Furthermore, Eq. (19) can be decomposed into a flavor nonsinglet (NS) and singlet (S) component

$$g_1^p(x,Q^2)=g_{1,NS}(x,Q^2)+g_{1,S}(x,Q^2)$$
(20)
where

$$g_{1,NS}(x,Q^2)=\frac{1}{2}\underset{q}{}(e_q^2e^2)(\delta q+\delta \overline{q})$$
(21)
with $`e^2=\frac{1}{f}_qe_q^2`$ (e.g. for f=3 light $`u,d,s`$ flavors $`e^2=\frac{2}{9}`$)and

$$g_{1,S}(x,Q^2)=\frac{1}{2}e^2\underset{q}{}(\delta q+\delta \overline{q})\frac{1}{2}e^2\delta \mathrm{\Sigma }.$$
(22)
where the sum usually runs over the light quark-flavors $`q=u,d,s`$, since the heavy quark contributions (c, b,…) have preferably to be calculated perturbatively from the intrinsic light quark ($`u,d,s`$) and gluon ($`g`$) partonic-constituent of the nucleon.

We are now in a position to present the results for the proton polarized structure function, $`g_1^p`$. In Eq. (19) we presented polarized proton structure function in a proton. In figure (4) the results are shown at some values of $`Q^2`$. The data points are from . For the neutron structure function one only needs to interchange $`u_v`$ and $`d_v`$ in the Eq. (19). For deuteron it is given by:

$$g_1^d(x,Q^2)=\frac{1}{2}[g_1^p(x,Q^2)+g_1^n(x,Q^2)](1\frac{3}{2}\omega _D)$$
(23)
where $`\omega _D0.058`$ accounts for the D-state admixture in the deuteron wave function. Quantities of particular significance are the first moment $`(n=1)`$ of polarized parton distributions $`\delta f(x,Q^2),`$

$$\mathrm{\Delta }f(Q^2)_0^1𝑑x\delta f(x,Q^2)\text{ , }f=q,\overline{q},g$$
(24)
We calculated first moments (total polarizations) $`\mathrm{\Delta }f(Q^2)`$ of the polarized parton distributions. These results shown in Table II.

According to Eqs. (17) and (18), $`\mathrm{\Delta }(q+\overline{q})`$ is the net number of right-handed quarks of flavor $`q`$ inside a right-handed proton and thus $`\frac{1}{2}\mathrm{\Delta }\mathrm{\Sigma }`$, is a measure of how much all quark flavors contribute to the spin of the proton. Similarly, $`\mathrm{\Delta }g`$ represents the total gluonic contribution to the spin of the nucleon.

Here we have defined the first moment of $`g_1`$ (Ellis-Jaffe sum rule) by

$$\mathrm{\Gamma }_1^{p,n,d}(Q^2)_0^1𝑑xg_1^{p,n,d}(x,Q^2)$$

Our values for $`\mathrm{\Gamma }_1^{p,n,d}`$ represents in Table III, which are in excellent agreement with the E143, E142, E154, SMC results from Ref. \[16-19\].

Figure (5) shows our results for $`xg_1^n`$ at $`Q^2=2,3,5`$ $`GeV^2`$ which are in agreement with the results of experimental data from .

In figure (6) we have shown $`xg_1^d`$ for some values of $`Q^2`$. The experimental data are from .

Finally by this model it is possible to attribute numbers $`\mathrm{}\mathrm{\Sigma }`$ and $`\mathrm{\Delta }g`$ to the quark and gluon spin content of the nucleons which describe their total (integrated) contribution to the nucleon spin in the following sense

$$\frac{1}{2}=\frac{1}{2}\mathrm{\Delta }\mathrm{\Sigma }(Q^2)+\mathrm{\Delta }g(Q^2)+L_z(Q^2)$$
(25)
where on the left hand side we have the spin $`(+\frac{1}{2})`$ of a polarized nucleon state and on the right hand side a decomposition in terms of $`\mathrm{\Delta }\mathrm{\Sigma }(Q^2)=\mathrm{\Delta }u+\mathrm{\Delta }d+\mathrm{\Delta }s+\mathrm{\Delta }\overline{u}+\mathrm{}`$, $`\mathrm{\Delta }g(Q^2)`$ and the relative orbital angular momentum $`L_z(Q^2)`$ among all the quarks and gluons. Unfortunately, the decomposition in Eq. (25) cannot directly be measured in experiments. Instead various other combinations of $`\mathrm{\Delta }\mathrm{\Sigma }`$ and $`\mathrm{\Delta }g`$ appear in experimental observables. Our calculation for study of nucleon spin presented in Fig. (7).

## 5 CONCLUSION

We have presented a model in which the polarized parton density is calculated. In this paper polarized valon distribution is derived from unpolarized valon distribution. After define the moment of polarized valon distribution and moments of all partons in proton, polarized parton density in a hadron are available. The results are used to evaluate the spin components of nucleon. It turns out that are results about polarized structure function agree well with all available experimental data on $`g_1`$ of proton ,neutron and deuteron. Also our prediction for $`\mathrm{\Gamma }_1`$ is good agreement with experimental data.

Table Caption

| $`v`$ | $`N`$ | $`\alpha `$ | $`\beta `$ | $`\gamma `$ | $`\eta `$ |
| --- | --- | --- | --- | --- | --- |
| $`U`$ | 0.0475 | -1.003 | -0.998 | 0.653 | -0.067 |
| $`D`$ | -0.015 | -1.247 | -0.951 | 8.351 | -2.063 |

Table I. The parameters defined in Eqs. (4,5) for $`U`$ and $`D`$ polarized valon distribution.

| $`Q^2(GeV^2)`$ | $`\mathrm{\Delta }u_v`$ | $`\mathrm{\Delta }d_v`$ | $`\mathrm{\Delta }\mathrm{\Sigma }`$ | $`\mathrm{\Delta }\overline{q}`$ | $`\mathrm{\Delta }g`$ |
| --- | --- | --- | --- | --- | --- |
| 2 | 0.609 | -0.306 | 0.329 | 0.004 | 0.876 |
| 3 | 0.609 | -0.306 | 0.334 | 0.005 | 0.992 |
| 5 | 0.609 | -0.306 | 0.341 | 0.006 | 1.139 |
| 10 | 0.609 | -0.306 | 0.350 | 0.008 | 1.334 |

Table II. The first moments of polarized parton distributions, $`\mathrm{\Delta }u_v`$, $`\mathrm{\Delta }d_v`$, $`\mathrm{\Delta }\mathrm{\Sigma }`$, $`\mathrm{\Delta }\overline{q}`$ and $`\mathrm{\Delta }g`$ for some value of $`Q^2`$.
$`Q^2(GeV^2)`$ $`\mathrm{\Gamma }_1^p`$ $`\mathrm{\Gamma }_1^n`$ $`\mathrm{\Gamma }_1^d`$    2 0.1212 -0.0313 0.0412  3 0.1219 -0.0307 0.0417  5 0.1225 -0.0299 0.0424  10 0.1236 -0.0289 0.0433
Table III. The first moments of proton, neutron, and deuteron polarized structure function for $`Q^2=2,3,5,10`$ $`GeV^2`$.


text:
# A dynamical mechanism for establishing apsidal resonance

## 1 Introduction

Orbital eccentricities in exo-planetary systems111http://www.exoplanets.org discovered thus far are often surprisingly large and have proven to be a major puzzle in understanding these systems. At least two systems with multiple planets in eccentric orbits ($`\upsilon `$ Andromedae and HD 83443) are suspected of exhibiting secular apsidal resonance (Chiang, Tabachnik & Tremaine 2001, Wu & Goldreich 2002). Apsidal resonance is the phenomenon of phase-locking of the apsidal longitudes of two orbits, such that the two planets have a common average rate of apsidal precession and the angular difference of their apsidal longitudes, $`\mathrm{\Delta }\varpi =\varpi _1\varpi _2`$, librates around 0.

We describe here a dynamical mechanism for establishing apsidal resonance in a pair of planets that are initially on nearly circular orbits. We use classical analysis from celestial mechanics to show that an impulse perturbation that imparts an eccentricity to one of the orbits, excites the other planet’s eccentricity on a secular timescale and also results in the libration of the relative apsidal longitude for a wide range of initial conditions. A plausible cause of an impulse perturbation is the ejection of a planet from the system; torques from an exterior primordial disk may also cause such an eccentricity impulse.

## 2 Secular dynamics

Consider a pair of planets in well-separated orbits, whose orbital period ratio is not too close to a ratio of two small integers. Under the assumption that the relative inclination of their orbital planes is small, the secular interactions in such a system preserve the orbital semimajor axes and induce quasi-periodic variations in the other orbital parameters. For small eccentricities and small relative inclinations, the eccentricity-apsidal line variations are decoupled from the inclination-nodal line variations (Murray & Dermott 1999). Under these conditions, the secular variations of the eccentricities and apsidal longitudes are described by a set of first order linear differential equations for the components of the so-called eccentricity vectors, $`(k_j(t),h_j(t))=e_j(\mathrm{cos}\varpi _j,\mathrm{sin}\varpi _j)`$, for each of the planets, where $`e_j`$ is the eccentricity and $`\varpi _j`$ is the apsidal longitude. The general solution for two planets is a superposition of two eigenmodes,

$`h_j(t)`$ $`=`$ $`E_j^{(1)}\mathrm{sin}(g_1t+\beta _1)+E_j^{(2)}\mathrm{sin}(g_2t+\beta _2),`$ (1)
$`k_j(t)`$ $`=`$ $`E_j^{(1)}\mathrm{cos}(g_1t+\beta _1)+E_j^{(2)}\mathrm{cos}(g_2t+\beta _2).`$ (2)
Let $`m_{}`$, $`m_1`$, $`m_2`$ be the masses of the star, the inner planet and the outer planet, respectively; $`n_1`$, $`n_2`$ the orbital mean motions of the inner and outer planet, respectively; $`\alpha =a_1/a_2<1`$ the ratio of the semimajor axes, $`b_s^{(l)}(\alpha )`$ a Laplace coefficient. We also define $`\mu =m_1/m_2`$, and $`B=b_{3/2}^{(2)}/b_{3/2}^{(1)}`$. With these definitions, the mode frequencies are

$$g_{1,2}=\frac{1}{8}\sqrt{\alpha }b_{3/2}^{(1)}(\alpha )\frac{m_2}{m_{}}n_2\left\{1+\mu \sqrt{\alpha }\pm \left[(1\mu \sqrt{\alpha })^2+4\mu \sqrt{\alpha }B^2\right]^{1/2}\right\},$$
(3)
and the eigenvector amplitudes $`𝐄^{(j)}`$ are given, up to a normalization factor, by

$$\frac{E_2^{(j)}}{E_1^{(j)}}=\frac{1\mu \sqrt{\alpha }[(1\mu \sqrt{\alpha })^2+4\mu \sqrt{\alpha }B^2]^{1/2}}{2B}\rho _j,$$
(4)
where the $`\mathrm{`}^{}`$ sign in $``$ corresponds to $`j=1`$, and the $`\mathrm{`}+^{}`$ sign to $`j=2`$. The phases $`\beta _j`$ and the normalization factors for the amplitudes of the modes are determined by initial conditions.

We assume initially nearly circular orbits, as is expected of planet formation in a dissipative protoplanetary disk. And we consider an impulse perturbation222In general, an impulse perturbation would result in not only an eccentricity impulse but also impulse perturbation to other orbital elements. The secular dynamics described here remains valid provided we use the post-impulse value of $`a_2`$ in evaluating $`\alpha `$ and $`n_2`$. that imparts a finite eccentricity, $`e_{2\mathrm{f}}`$, to the outer planet’s orbit, at time $`t=0`$. (The case of an impulse to the inner planet’s orbit is solved in a similar way.) Then, with initial conditions $`h_1(0)=h_{10},k_1(0)=k_{10},h_2(0)=0,k_2(0)=e_{2\mathrm{f}}`$, we have

$`E_1^{(1)}\mathrm{cos}\beta _1`$ $`=`$ $`{\displaystyle \frac{\rho _2k_{10}e_{2\mathrm{f}}}{\rho _2\rho _1}},E_1^{(1)}\mathrm{sin}\beta _1={\displaystyle \frac{\rho _2h_{10}}{\rho _2\rho _1}},`$ (5)
$`E_1^{(2)}\mathrm{cos}\beta _2`$ $`=`$ $`{\displaystyle \frac{\rho _1k_{10}e_{2\mathrm{f}}}{\rho _1\rho _2}},E_1^{(2)}\mathrm{sin}\beta _2={\displaystyle \frac{\rho _1h_{10}}{\rho _1\rho _2}}.`$ (6)
This solution admits either circulation or libration of the relative apsidal longitude, $`\mathrm{\Delta }\varpi `$. Examples of three types of possible behavior are shown in Fig. 1. In each of these cases, the cumulative probability density function of $`|\mathrm{\Delta }\varpi |`$ and of $`e_1`$ (i.e., the fraction of time that $`|\mathrm{\Delta }\varpi |,e_1`$ is less than some value) is also shown.

In the case of zero initial eccentricity, $`e_{10}=0`$ \[Fig. 1(a–d)\], mode 1 and mode 2 appear with equal amplitude but opposite phase in the inner planet’s eccentricity vector. The time evolution of the planets’ eccentricities and of $`\mathrm{\Delta }\varpi `$ is given by

$`e_1(t)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}e_{2\mathrm{f}}}{|\rho _1\rho _2|}}\left[1\mathrm{cos}(g_1g_2)t\right]^{1/2},`$ (7)
$`e_2(t)`$ $`=`$ $`{\displaystyle \frac{e_{2\mathrm{f}}}{|\rho _1\rho _2|}}\left[\rho _1^2+\rho _2^22\rho _1\rho _2\mathrm{cos}(g_1g_2)t\right]^{1/2},`$ (8)
$`\mathrm{tan}\mathrm{\Delta }\varpi (t)`$ $`=`$ $`{\displaystyle \frac{\rho _2\rho _1}{\rho _2+\rho _1}}\mathrm{tan}[{\displaystyle \frac{1}{2}}(g_1g_2)t+{\displaystyle \frac{\pi }{2}}].`$ (9)
The eccentricity of the inner planet grows from its initial value of zero to a value comparable to $`e_{2\mathrm{f}}`$ on the secular timescale $`|g_1g_2|^1`$; the eccentricities vary periodically (but not sinusoidally) with a common frequency, $`g_1g_2`$, and with opposite phases. The inner planet’s eccentricity has a maximum excursion to $`2e_{2\mathrm{f}}/|\rho _1\rho _2|`$, while the eccentricity of the outer planet varies between a minimum of $`e_{2\mathrm{f}}|\rho _1+\rho _2|/|\rho _1\rho _2|`$ and a maximum of $`e_{2\mathrm{f}}`$. The apsidal longitude difference is limited to $`(90^{},90^{})`$; it oscillates with the same frequency, $`g_1g_2`$. The singularity of the tangent function at $`\pm \frac{1}{2}\pi `$ implies a discontinuity in the time evolution of $`\mathrm{\Delta }\varpi `$. This discontinuity is not unphysical (see Fig. 2 below): as the eccentricity vector of the inner planet passes periodically through the origin, its apsidal longitude changes discontinuously.

A small non-zero initial eccentricity, $`e_{10}e_{2\mathrm{f}}`$, results in apsidal libration with an amplitude smaller than 90 for values of $`\varpi _{10}`$ that are approximately within $`\pm 90^{}`$ of $`\varpi _{20}`$ \[e.g., Fig. 1(e–h)\], and apsidal circulation otherwise \[e.g., Fig. 1(i–l)\]. In all cases, the secular dynamics following an impulse perturbation to $`e_2`$ yields a low probability of finding $`|\mathrm{\Delta }\varpi |>90^{}`$ \[see Fig. 1(c,g,k)\], and only $`1/3`$ probability of finding $`e_1`$ less than half its maximum \[see Fig. 1(d,h,l)\]. As we show below, a small non-zero $`e_{10}`$ results in apsidal resonance with approximately 50% probability.

It is illustrative to consider the limiting case $`m_10`$. In this case, the orbital elements of the outer planet are constant: $`k_2(t)=e_{2\mathrm{f}},h_2(t)=0`$; for the inner planet, we have the following solution:

$`k_1(t)`$ $`=`$ $`+(k_{10}e_{1\mathrm{f}})\mathrm{cos}gt+h_{10}\mathrm{sin}gt+e_{1\mathrm{f}},`$ (10)
$`h_1(t)`$ $`=`$ $`(k_{10}e_{1\mathrm{f}})\mathrm{sin}gt+h_{10}\mathrm{cos}gt,`$ (11)
with

$$e_{1\mathrm{f}}=Be_{2\mathrm{f}}\text{and}g=\frac{1}{4}n_1\frac{m_2}{m_{}}\alpha ^2b_{3/2}^{(1)}(\alpha ).$$
(12)
This solution admits a simple geometrical interpretation as illustrated in Fig. 2a: the inner planet’s eccentricity is the vector sum of a ‘forced’ component of magnitude $`e_{1\mathrm{f}}`$ which is aligned with the outer planet’s apsidal longitude (taken to be zero here, without loss of generality), and a ‘free’ component of magnitude

$$e_{1,\mathrm{free}}=(e_{10}^2+e_{1\mathrm{f}}^22e_{1\mathrm{f}}k_{10})^{1/2}$$
(13)
which rotates with frequency $`g`$. Note that $`e_{1\mathrm{f}}`$ is independent of the mass of the outer planet, but the frequency $`g`$ is proportional to that mass.

Apsidal libration occurs if $`e_{1,\mathrm{free}}e_{1\mathrm{f}}`$. In the $`(k_{10},h_{10})`$ plane, this condition is satisfied in the area of intersection of two circles, one of radius $`e_{10}`$ centered at the origin, the other of radius $`e_{1\mathrm{f}}`$ centered at $`(e_{1\mathrm{f}},0)`$, as shown in Fig. 2b. Assuming random values of $`\varpi _{10}`$ and $`e_{10}\epsilon `$, the probability $`P_{\mathrm{lib}}`$ that the impulse on the outer planet will result in apsidal libration is simply the intersection area as a fraction of the area of the circle of radius $`\epsilon `$. This is readily evaluated:

$$P_{\mathrm{lib}}(e_{10}\epsilon )=\frac{1}{\pi }\left[\mathrm{cos}^1\left(\frac{\epsilon }{2e_{1\mathrm{f}}}\right)+2\left(\frac{e_{1\mathrm{f}}}{\epsilon }\right)^2\mathrm{sin}^1\left(\frac{\epsilon }{2e_{1\mathrm{f}}}\right)\frac{e_{1\mathrm{f}}}{\epsilon }\sqrt{1\left(\frac{\epsilon }{2e_{1\mathrm{f}}}\right)^2}\right].$$
(14)
For $`\epsilon e_{2\mathrm{f}}`$, we have $`P_{\mathrm{lib}}=0.5`$: nearly all values of $`\varpi _{10}`$ between $`90^{}`$ and $`+90^{}`$ result in apsidal resonance. For $`\epsilon =e_{1\mathrm{f}}`$, the probability of apsidal resonance drops only slightly, to $`0.4`$.

For the general case, $`m_10`$, the probability of apsidal resonance can be determined numerically. It is easy to see from the general solution (eq. 4) that $`P_{\mathrm{lib}}`$ has only a weak dependence on the mass ratio $`\mu =m_1/m_2`$, and $`P_{\mathrm{lib}}`$ remains close to 0.5.

When apsidal libration occurs, the libration amplitude is a steep function of $`\varpi _{10}`$ and $`e_{10}/e_{2\mathrm{f}}`$. This is illustrated in Fig. 3, for one choice of parameters, $`\mu =0.5`$ and $`\alpha =0.33`$, corresponding to the case of planets C and D of $`\upsilon `$ Andromedae.

## 3 Discussion

We have shown that in a two planet system, an impulse perturbation on the eccentricity of the outer planet can excite the eccentricity of the inner planet on a secular timescale, and will result in apsidal libration with a $`50\%`$ probability. There are two significant characteristics of the secular dynamics following an eccentricity impulse, whether or not it results in apsidal resonance: (i) the inner planet’s eccentricity has large amplitude variation, and it drops to its initial small value periodically; (ii) only a small fraction of time is spent at large values of $`|\mathrm{\Delta }\varpi |`$ and small values of $`e_1`$ (see Fig. 1). Therefore, observing small values of $`|\mathrm{\Delta }\varpi |`$ and large values of $`e_1`$ is not surprising if the system has suffered an impulse perturbation in its history and the secular dynamics indicates large amplitude variations.

We note that in the opposite limit of a slow adiabatic perturbation which increases one planet’s eccentricity on a timescale much longer than the secular timescale, apsidal resonance will occur with nearly 100% probability and the resulting libration amplitude will be very small (for initially circular orbits). It is rather remarkable that the secular dynamics of two planets admits a high probability of apsidal resonance in both the impulse and the adiabatic limit of eccentricity perturbation of one planet. We defer detailed analysis to a future study, but we note here that the apsidal libration amplitude provides an observational diagnostic: large amplitudes favor the impulse mechanism whereas small amplitudes favor the adiabatic mechanism.

One possible mechanism for the eccentricity impulse is planet-planet scattering. Highly eccentric orbits can be produced by gravitational interactions of planets initially on circular orbits that are close to the threshold of orbital instability. This has been proposed as a possible explanation for the eccentric orbits of exo-planets (Rasio & Ford 1996, Weidenschilling & Marzari 1996, Lin & Ida 1997, Ford et al. 2001, Marzari & Weidenschilling 2002). There are three possible outcomes of such interactions: a final stable configuration different from initial, the merger of two planets, or the ejection of one planet. Numerical experiments thus far have explored only a small fraction of the parameter space of two and three giant-planet systems. In systems of two equal mass planets, comparable branching ratios for the three possible outcomes are found; the two former outcomes do not result in eccentricity excitation, but the third – ejection of a planet – does, with typical eccentricity of the surviving planet in the range 0.4–0.8 (Ford et al. 2001). One problematic aspect is that two planet systems must have initial orbits fine-tuned close to the threshold of dynamical instability, i.e. orbital separation close to $`2\sqrt{3}`$ times their mutual Hill radius (Gladman 1993). But this problem does not exist in systems of three or more planets, where chaotic instabilities can occur over a wider range of orbital separations; such systems can persist in quasi-stable orbital configurations for long periods of time, possibly exceeding $`10^8`$ years, before becoming dynamically unstable; then the most common outcome is the ejection of one planet (Marzari & Weidenschilling 2002). Also, in the more realistic case of unequal masses, the probability of planet ejection is likely to be enhanced further.

We consider the ejection into an unbound orbit of a third planet by gravitational scattering from the outer of the two planets whose secular dynamics we have analyzed in section 2. (For comparable masses, ejection by scattering from the outer of the two planets is more likely than from the inner planet because the escape velocity is smaller at larger distance from the star: $`v_{\mathrm{esc}}r^{1/2}`$; it is also likely to produce the least direct perturbation to the inner members of the system.) Assuming that the ejected planet leaves on a nearly parabolic orbit, energy and angular momentum conservation yield the following estimate for the final eccentricity of the bound planet:

$$1e_{2\mathrm{f}}^2(1+\frac{m_3}{m_2}\frac{a_2}{a_3})\left[1+\frac{m_3}{m_2}\sqrt{\frac{m_{}+m_2}{m_{}+m_3}}\left(\sqrt{\frac{a_3}{a_2}}\sqrt{\frac{2q_3}{a_2}}\right)\right]^2,$$
(15)
where $`m_3,a_3,q_3`$ are the mass, initial semimajor axis and final periastron distance of the ejected planet; $`a_2`$ is the initial semimajor axis of $`m_2`$, and is related to its final semimajor axis by energy conservation: $`a_{2\mathrm{f}}m_2a_2a_3/(m_2a_3+m_3a_2)`$. Unfortunately, there is no additional dynamical constraint or simple argument to constrain the final periastron distance of the ejected planet. So the final eccentricity of the bound planet cannot be predicted analytically and must be determined numerically. Still, we can make the following estimate: since the periastron distance of the ejected planet is likely to be not greatly in excess of the apoastron distance of the surviving planet, then, for $`q_3/a_2`$ 1–1.5, the final eccentricity of the bound planet will be in the range of 0.2–0.7 for $`m_3/m_2`$ 0.1–0.7.

Once a dynamical instability sets in, the subsequent evolution is highly chaotic and unpredictable in detail. A simplified description, based upon the presently available numerical simulations, is as follows: at first, close encounters between $`m_2`$ and $`m_3`$ lead to a rapid increase in their orbital eccentricities, over $`1010^2`$ orbital periods; this is followed by a longer period, $`10^4`$ orbits, during which many weak distant encounters gradually and stochastically increase the orbital period and apoastron of $`m_3`$, the planet-to-be-ejected. In the present context, the initial eccentricity excitation of the surviving planet, $`m_2`$, satisfies the impulse approximation as the timescale of the first stage is much shorter than the secular timescale (which is of order $`(m_{1,2}/m_{})^1`$ times the orbital period, or $`10^3`$ orbital periods for jovian mass planets). Thus, the secular dynamics of the two surviving planets, $`m_1`$ and $`m_2`$, would be as described in Section 2. The weak distant encounters that eventually lead to the ejection of the third planet would produce small perturbations to the secular solution, as would the gravitational interactions between $`m_1`$ and $`m_3`$. (Neglecting the latter is justified, in the lowest order, as the $`m_1`$,$`m_3`$ separation is much greater than the $`m_2`$,$`m_3`$ separation.) The effects of these perturbations on the secular dynamics presented above will be considered in a future study.

The case of $`\upsilon `$ Andromedae The analysis given in the previous section can be used to estimate the pre-impulse orbital eccentricity of at least one of the planets. A good candidate amongst the known exo-planetary systems to apply this theory is the $`\upsilon `$ Andromedae system, which represents perhaps the most (dynamically) constrained of the known exo-planetary systems. A recent orbital solution for this system (as quoted in Chiang et al. 2001) gives present eccentricities of planets C and D of 0.25 and 0.34, respectively, and yields secular dynamics for these two planets similar to that shown in Fig. 1(e,f), with apsidal libration amplitude of 25–35 (for relative orbital inclinations $`20^{}`$). An independent orbital solution obtained previously by Stepinski et al. 2000 (from an older and therefore smaller set of observational data) has orbital parameters within $`2\sigma `$ uncertainties of the more recent solution; it yields secular dynamics similar to that shown in Fig. 1(a,b), with apsidal libration amplitude of 80–90. The libration period for this system is $`7\times 10^3`$ yr. Now, if we allow that planet D suffered an eccentricity impulse, $`e_D0.35`$, in its history, we can estimate that the initial eccentricity of planet C could have been as small as $`0.07`$ if the apsidal libration amplitude is 25–35 deg, but even smaller, 0–0.01, if the amplitude is 80–90.

After submitting this paper, we became aware of a recent preprint by Chiang & Murray (2002) who propose an adiabatic eccentricity perturbation to explain the apsidal resonance in the $`\upsilon `$ Andromedae planetary system. They invoke torques from an exterior massive primordial disk to provide the adiabatic eccentricity excitation to the outermost planet D. However the adiabaticity of the perturbation is owed to a specific choice of disk parameters, which are not well constrained; a different, and perhaps equally plausible, choice of disk parameters may instead provide an impulse perturbation (E. Chiang 2002, personal communication); then the secular dynamics described here would apply.

As mentioned previously, the apsidal libration amplitude provides an observational diagnostic: large amplitudes favor the impulse mechanism whereas small amplitudes favor the adiabatic mechanism. In the case of $`\upsilon `$ Andromedae, there remain significant uncertainties in the orbital parameters, and at present it appears difficult to discriminate between the adiabatic and the impulse mechanisms. The 25–35 deg apsidal libration amplitude found by Chiang et al. (2001) is neither clearly small nor clearly large; in both the adiabatic mechanism and the impulse mechanism, an initial eccentricity near 0.06–0.07 of planet C can account for this libration amplitude. On the other hand, the large libration amplitude of 80–90 deg found by Stepinski et al. (2000) would be explained by initial $`e_C0.01`$ in the impulse limit, but cannot be explained naturally in the adiabatic limit. We urge further observations and analysis to improve the accuracy and fidelity of orbital solutions of exo-planetary systems; this would help to constrain their dynamical history.

We acknowledge comments from E. Chiang, discussions with F. Rasio, and research support from NASA.


text:
# 1 Introduction

## 1 Introduction

Atomic clusters are aggregates of atoms held together by the same forces that cause, for example, phase transition from vapor to liquid, formations of crystals, etc. Cluster sizes range from as few as three atoms up to several hundred atoms. The physical and chemical characteristics of a cluster often varies with its size. In fact, even the addition of a single atom can result in an entirely different structure. Only by successively adding more and more atoms will a crystal-like structure eventually be produced and some knowledge of the condensed phase attributes be determined .

The study of atomic clusters has steadily been increasing over the past decade . Of particular interest is the cluster conformation (structure) which has the lowest total internal energy $`E`$. Knowledge of this minimum energy conformation provides valuable clues relating to the chemical and physical properties of the cluster. Unfortunately, searching for the globally minimum energy state of a cluster has proven to be enormously difficult. Indeed, Wille and Vennik  showed that locating the globally minimum energy state of a cluster of identical atoms—the homonuclear case—belongs in the class of NP-hard problems. This means there is little hope of exactly solving the problem in finite time for even moderate cluster sizes.

The purpose of this paper is two-fold. First, it will be shown why existing homonuclear proofs, and work from other related problems, cannot be used for the heteronuclear case where not all of the atoms are identical. Secondly, a proof will be presented which does show solving the heteronuclear problem is NP-hard.

## 2 Preliminaries

The problem of finding this globally minimum energy structure is equivalent to optimizing $`E`$ with respect to variations in all $`3N6`$ degrees of freedom where $`N`$ is the cluster size. One method of solving this optimization problem is to explore the potential energy surface (PES) composed of all possible cluster conformations. Each point on this surface represents a unique spatial arrangement of the constituent atoms. This multidimensional surface is characterized by numerous local minima, each indicating an energetically stable structure. If it were possible to enumerate all of these minima—and the saddles that link them—we could, in principle, describe the dynamics of chemical reactions governed by that surface. Unfortunately, enumerating all minima is extremely difficult because of their large number. Berry  indicates the number of geometrically distinct minima tends to grow exponentially with $`N`$. Moreover, the number of permutational isomers grows factorially with $`N`$.

A number of general methods have been proposed for finding global minima on a PES, in particular, and on hypersurfaces in general. For example, Monte Carlo methods , eigenvector following , evolution computation techniques -, lattice optimization/relaxation techniques , and PES deformation techniques  have all been used with varying degrees of success. After formally defining the clustering problem, we will discuss some of these techniques in greater detail.

The potential energy function for a cluster of $`N`$ atoms is given by

$$V(r^N)=\frac{1}{2}\underset{\begin{array}{c}i,j\\ ij\end{array}}{\overset{N}{}}v(r_ir_j)$$
(1)
where $`r^N=(r_1,r_2,\mathrm{},r_N)`$, $`r_i`$ is the position vector of the $`i`$th atom and $`v(r_ir_j)`$ is a function representing the pairwise interaction between atoms. Lennard-Jones or Morse functions are frequently used for these functions. Our optimization problem of interest is the Discrete Cluster Problem (DCP) which is formally defined as follows:111The definition given covers both homonuclear and heteronuclear systems.

DISCRETE CLUSTER PROBLEM

INSTANCE: Given finite number of points in real space, a distance $`d(i,j)Z^+`$ between two points $`i,j`$, an integer $`NZ^+`$ and a known potential energy function defined by (1).

QUESTION: Is there a way to assign $`N`$ atoms to $`N`$ unique points so as to minimize the sum of their pairwise interactions?

The potential energy of a cluster actually equals the sum of $`N`$-body interactions. Restricting this sum to only two-body interactions provides only a *qualitative* approximation. Three-body interactions can be important, but are usually ignored in first order approximations . Although many-body interactions are needed for *quantitative* modeling, little is actually known about higher order terms. Many-body potential energy functions do exist , though in practice, only two-body terms are used for the sake of computational speed. Consequently, the discussion here will likewise be restricted to the two-body case.

It may appear that the large amount of work done with hard-sphere packing problems will be helpful in solving instances of DCP. Unfortunately, such is not the case because the objective of the two problems are quite different. Hard-sphere packings try to place $`N`$ spheres in Euclidean space so that all can, without overlap, fit within as small a volume as possible . DCP deals with soft, compliant spheres which interact via pairwise interaction functions. Figure 1 shows a Lennard-Jones function, which is typical. Therein lies the major difference between hard-sphere and soft-sphere systems: the former has no preferred distance between the spheres while the latter does. Put another way, a hard-sphere packing algorithm attempts to minimize the interatomic distance $`r`$. Yet, a comparison with Figure 1 clearly shows this does not yield the lowest energy state for an atomic pair. Hard-sphere packing studies can thus be expected to provide little help.

An algorithm that searches for solutions to the homonuclear version of DCP was recently proposed by Hendrickson . The general approach exploits structural information to decompose the global optimization problem into a set of smaller optimization problems. More specifically, $`N`$ objects in Euclidean space are represented as an undirected graph where the edge weights reflect the interobject distance. A subgraph is identified, the relative positions of the vertices are optimized, and the subgraph is then treated as a rigid body. These rigid bodies are ultimately combined to determine the overall structure. Such an approach won’t work for the heteronuclear case, but this discussion will be differed until Section 3.

Northby  performed a lattice based search that takes an initial conformation and allows it to relax to an energy minimum. Essentially, a random search is conducted to compile a list of isomers. After pruning any geometrically equivalent isomers, the remaining conformations are relaxed under a Lennard-Jones pairwise interaction function. While this technique will work with heteronuclear clusters, it does depend upon a gradient search which makes it susceptible to stopping at local optima.

Kostrowicki et al. , presented a technique that “deforms” the PES in such a way that the number of local minima is dramatically reduced. This deformation requires solving a partial differential equation called a diffusion equation. All atoms are assumed to interact by a Lennard-Jones type of function which is approximated by a linear combination of Gaussian functions to permit expressing the answer analytically. Solving these diffusion equations is not necessarily easy as there are some pathological conditions. For example, numerical problems can occur when the atoms are separated by large distances and so interact only weakly. This technique works best if one can suitably modify the boundary conditions for the diffusion equation, but care must be taken to ensure the minimum energy state is not removed from the deformed PES.

Wales and Doye  described another PES transformation technique called “basin hopping”. Here the PES is converted into a set of basins of attraction for the local minima. This approach does has the advantage of not altering the energies of the minima. It has been used to find many of the local minima for Lennard-Jones cluster sizes up to $`N=110`$, The authors suggest some improvements that could be made to improve the efficiency of their approach.

*It is important to note that none of these techniques will, with certainty, guarantee a successful search for the globally minimum energy configuration*. In the next section we will see why solving instances of DCP has proven to be so difficult.

## 3 Complexity in Homonuclear Clusters

An instance of DCP can be solved by exploring the PES associated the cluster, and any algorithm that does this can be used with both homonuclear and heteronuclear systems. A number of researchers indicate the number of local optima in a PES grows exponentially with $`N`$ . So, exactly how much effort is required to explore an exponentially large hypersurface? Some insight can be found from the theory of NP-completeness .

Suppose each time step a new point on the hypersurface is visited. Repeating this process until all points are visited—a procedure guaranteed to find the minimum energy state—will take an exponential number of time steps. Now consider an arbitrary NP-complete problem $`𝒫`$. If we found an algorithm that could solve $`𝒫`$ in polynomial time, then this algorithm could solve every NP-complete problem in polynomial time. Unfortunately, no algorithm that solves an NP-complete problem has ever been found which executes in less than an exponential number of steps. This suggests the effort required to conduct a blind search of the hypersurface is similar to the effort required to solve an NP-complete problem. NP-complete problems are computationally intractable, which gives some idea of the difficulty one faces in solving DCP. (NOTE: This does *not* prove DCP is NP-complete. It merely shows the ramifications of having so many local optima.)

NP-complete problems are defined as decision problems, that is, problems answered by either ‘yes’ or ‘no’. NP-hard problems ask for the optimal solution to an NP-complete problem. And, they have at least the same level of difficulty to solve as does the corresponding NP-complete problem. One particular NP-hard problem plays a pivotal role in the complexity analysis of DCP. It is the well known Traveling Salesman Problem (TSP) which is formally defined as follows:

TRAVELING SALESMAN PROBLEM

INSTANCE: A finite set $`C=\{c_1,c_2,\mathrm{},c_m\}`$ of cities, and a distance $`d(c_i,c_j)Z^+`$ for each pair of cities $`c_i,c_jC`$.

QUESTION: What permutation

$$[c_{\pi (1)},c_{\pi (2)},\mathrm{},c_{\pi (m)}]$$

of $`C`$ will minimize the tour length

$$\left\{\underset{i=1}{\overset{m1}{}}d(c_{\pi (i)},c_{\pi (i+1)})\right\}+d(c_{\pi (m)},c_{\pi (1)})\mathrm{?}$$

Wille and Vennik  used TSP to prove solving a homonuclear instance of DCP is NP-hard. It is worthwhile to examine their proof in some detail since our proof of NP-hardness for the heteronuclear case follows a similar line of reasoning.

Instances of DCP can be expressed in graph-theoretical terms using a graph $`G=(V,E)`$ with vertex set $`V`$ and edge set $`E`$. The graph $`G`$ is complete (i.e., the edge $`e=u,vEu,vV`$). Furthermore, $`|V|N`$. The edges are assigned weights $`w(e)eE`$ where the weights reflect the interaction between vertices. An instance of DCP is therefore equivalent to selecting $`V^{}V`$ with $`|V^{}|=N`$ so that

$$\frac{1}{2}\underset{\begin{array}{c}u,vV^{}\\ e=u,v\end{array}}{}w(e)=\text{minimal}$$
(2)
The proof of Wille and Vennik uses an undirected graph $`G=(V,E)`$ with $`|V|=N(N1)`$ vertices (see Figure 2). Each vertex is labeled with ($`c_i,c_j`$) where $`i,j=1,\mathrm{},N`$ and $`ij`$. This indicates city $`c_j`$ is visited immediately after city $`c_i`$. Edges are then unordered pairs of the form

$$e=(c_i,c_j)(c_k,c_l)$$

with edge weights

$$w(e)=\{\begin{array}{cc}0\hfill & \text{if}ik,jl\hfill \\ +\mathrm{}\hfill & \{\begin{array}{c}\text{if}i=k,jl\hfill \\ \text{if}j=l,ik\hfill \\ \text{if}l=i,k=j\hfill \end{array}\hfill \\ d(c_k,c_l)\hfill & \text{if}l=i,kj\hfill \\ d(c_i,c_j)\hfill & \text{if}k=j,li\hfill \end{array}$$
(3)
Selecting $`N`$ vertices, $`\{(c_{\pi (i)},c_{\pi (i+1)});i=1,\mathrm{},N,c_{\pi (N+1)}=c_{\pi (1)}\}`$, so that the sum of the weights is minimal is equivalent to finding a minimal length tour $`c_{\pi (1)},c_{\pi (2)},\mathrm{},c_{\pi (N)}`$ thus solving an instance of TSP. By restriction , this also makes DCP NP-hard to solve. This completes the Wille and Vennik proof.

The weight assignments given in Eq. (3) require some clarification. In TSP, weight is equivalent to distance whereas in DCP weight is equivalent to pairwise potential energy. Any plausible solution to an instance of TSP can visit each city but one time. In valid tour moves, the edge weight equals the intercity distance. Disjointed tours have edges with zero weight—effectively removing that edge from the edge set $`E`$. This insures all tours can be defined as a permutation of $`m`$ cities (see below). An illegal tour move has an edge with infinite weight. The total weight of the edges traversed in a tour measures the “goodness” of that tour; good tours have lower total edge weights. More specifically,

1. $`ik`$ and $`jl`$. This defines a disjointed tour where the tour visits $`c_i`$ followed by $`c_j`$ and $`c_k`$ followed by $`c_l`$. However, it is not shown what other cities were visited between $`c_j`$ and $`c_k`$. Therefore, it is not possible to describe the permutation and compute its tour length.
2. $`i=k`$ and $`jl`$. This defines an illegal tour that leaves the same city to visit two different cities.
3. $`j=l`$ and $`ik`$. This defines an illegal tour where two distinct cities visit the same next city.
4. $`l=i`$ and $`k=j`$. This defines an illegal cyclic tour among only two cities.
5. $`l=i,kj`$ or $`k=j,li`$. These are legal tours.

Each of the $`N`$ vertices in the minimal length tour is assigned a distinct atom from the cluster. In homonuclear clusters *all* one-to-one assignments of atoms to these $`N`$ vertices are equivalent since the edge weights are based only on interatomic distance. This begins to explain why the Wille and Vennik proof  and the Hendrickson algorithm  cannot be extended to the heteronuclear case.

First, consider a homonuclear cluster with $`N>2`$ atoms and suppose a new cluster is formed by swapping the spatial position of two atoms. It is not possible to tell any physical difference between the two clusters because all atoms are identical and the interatomic distances remain unchanged. Indeed, the new cluster will have a total energy identical to the original cluster because the pairwise interaction functions ($`v(r_ir_j)`$ in Eq. (1)) remain unchanged. Now suppose the cluster is composed of two distinct atom types, say Ar and Xe. The repulsive and attractive forces experienced by an Ar-Ar atom pair differs from those experienced by Ar-Xe or Xe-Xe pairs *even if all pairs are separated by the same interatomic distance* . Swapping atom positions in heteronuclear clusters changes the type of atoms which interact, altering the individual interaction functions, and giving a different total energy to the new cluster.

With homonuclear clusters it is acceptable to consider atoms as simple identical spheres where only interatomic distances contribute to the total energy. It is natural to model this system as an undirected graph where the edge weights reflect forces derived solely from the interatomic distances. The search algorithm from  takes this approach thereby permitting the cluster to be optimized in stages by optimizing the relative positions in subgraphs. The complexity proof given in  also made that assumption. In fact, the graph used for that proof was constructed specifically without requiring any pair type information to set the edge weights. That restriction was necessary to establish an equivalence between TSP and DCP.

In heteronuclear systems, both atom type and distance determine pairwise forces so the corresponding graph must have edge weights that take both distance and atom type into consideration. Even the relative positions of vertices from a subgraph cannot be optimized without the weights being set in this manner. Consequently, search algorithms such as  cannot be used for heteronuclear clusters because the overall cluster structure depends on connecting rigid bodies, formed from optimized subgraphs, where the edge weights only consider interatomic distance.

To apply the proof of Wille Vennik  to heteronuclear clusters, the graph would have to be augmented with additional vertices. For example, suppose in the original (homonuclear) graph two vertices $`i`$ and $`j`$ are connected by an edge $`e=i,j`$. Then $`w(e)`$ is based solely on the interatomic distance between atoms $`i`$ and $`j`$. The augmented graph would have a new vertex $`j^{}`$ where the added edge $`e^{}=i,j^{}`$ has a weight $`w(e^{})`$ computed from the interaction of two dissimilar atoms, spaced at a distance $`d(i,j^{})=d(i,j)`$. However, now the mere selection of $`N`$ vertices—sufficient to solve a homonuclear DCP—does not guarantee the correct mixture of atom types present in the heteronuclear cluster. This restricts the proof from  to only homonuclear systems.

## 4 Complexity in Heteronuclear Clusters

A different proof of complexity is needed for the heteronuclear clusters. This new proof makes use of the following NP-hard problem :

TRAVELING SALESMAN EXTENSION (TSE)

INSTANCE: A finite set $`C=\{c_1,c_2,\mathrm{},c_m\}`$ of cities, a distance $`d(c_i,c_j)Z^+`$ for each pair of cities $`c_i,c_jC`$, a bound $`BZ^+`$, and a “partial” tour

$$\mathrm{\Theta }=c_{\pi (1)},c_{\pi (2)},\mathrm{},c_{\pi (K)}$$

of $`K`$ distinct cities from $`C`$, $`1Km`$.

QUESTION: Can $`\mathrm{\Theta }`$ be extended to a full tour

$$c_{\pi (1)},c_{\pi (2)},\mathrm{},c_{\pi (K)},c_{\pi (K+1)},\mathrm{},c_{\pi (m)}$$

having total length $`B`$ or less?

Consider a heteronuclear cluster which has a single atom of one type (denoted by $`\alpha `$) and $`N1`$ atoms of a different type (denoted by $`\beta `$). A completely connected graph is $`G=(V,E)`$ with $`|V|=N(N1)`$ is constructed with each vertex labeled as described in Section 3. Without loss in generality preassign an $`\alpha `$-atom to all vertices with labels ($`c_1,c_j`$), $`j=2,\mathrm{},N`$ and preassign a $`\beta `$-atom to all remaining vertices.

The search begins with a scan of all edges that touch vertices with $`\alpha `$-atom assignments. Select the minimum weight edge.222As before, edges with weight 0 are effectively removed from the graph. This edge defines a partial tour $`c_{\pi (1)},c_{\pi (2)}`$. Now select $`N2`$ more vertices, $`\{(c_{\pi (i)},c_{\pi (i+1)});i=3,\mathrm{},N,c_{\pi (N+1)}=c_{\pi (1)}\}`$, so that the sum of these weights is minimal. (This search is limited to vertices with $`\beta `$-atom assignments in order to maintain the proper mixture of atom types in the corresponding cluster.) Finding a minimal length full tour is thus equivalent to solving an instance of TSE. This proves that DCP is NP-hard to solve for heteronuclear clusters as well.

## 5 Discussion

The complexity in solving DCP forces researchers to use heuristic search algorithms. Hill-climbing algorithms are not expected to do well because the PES has many local optima and it is highly likely the algorithm will quickly stop at one of them. Conversely, stochastic search algorithms can be quite effective in such multimodal hypersurfaces. It is therefore natural to ask if any one particular algorithm stands out as doing especially well against DCP.

This question is not easy to answer and we need to turn to the No Free Lunch (NFL) Theorem  for an explanation. Essentially, this theorem says if some optimization (search) algorithm performs particularly well over a certain class of optimization problems, then it most likely will not perform as well over all remaining optimization problems. This means one cannot choose, for example, simulated annealing to use for DCP just because it happens to work well for, say, scheduling problems. Direct comparisons between algorithms are insightful. But, without a conscientious attempt to make the comparisons fair, the results may be inconclusive, or in the worst case, be completely wrong . Monte Carlo techniques have been dominant in the area of cluster studies, which is why newly proposed search algorithms are normally compared against them.

Most notable among the new algorithms are the evolutionary algorithms, which conduct searches that mimic Darwinian evolution: a “population” of clusters evolve to a low energy state by altering existing configurations via stochastic reproduction operators. Natural selection determines which configurations survive to undergo further reproduction operations.

Comparisons between evolutionary algorithms and Monte Carlo techniques have favored the former, although such comparisons sometimes lack sufficient mathematical rigor. For example, Zeiri  concluded that a genetic algorithm converges faster than simulated annealing after comparing the average of five runs—an unusually small sample size. Normally results should be averaged over a considerably higher number of runs to help remove any potential bias in the random number generators used in the algorithms. Nevertheless, there is growing empirical evidence that says evolutionary algorithms consistently outperform the Monte Carlo techniques when applied against DCP  -. It appears as though evolutionary algorithms would be a good first choice search algorithm for cluster studies.

## Acknowledgement

The author wishes to thank the anonymous reviewers who made several pertinent and valuable suggestions.


text:
# Magnetic behavior of a non-extensive S-spin system: possible connections to manganites

## I Introduction

In 1988, TsallisTsallis (1988) proposed a $`q`$-dependent entropy functional which generalized the standard Maxwell-Boltzmann definition, to include non-extensive systems. In this formalism, "$`q`$" is called the *entropic index*, and measures the degree of nonextensivity of the system. A few years latter, Curado and TsallisCurado and Tsallis (1991) revised the formalism introducing the *unormalized* constraint to the internal energy and showed that from this entropy functional a nonextensive thermodynamics could be derived. This generalization included the classical thermodynamics (Maxwell-Boltzmann), which could be recovered when $`q`$ is set to unity.

After Curado and Tsallis work, this formalism has been successfully applied to various physical systems, where Maxwell-Boltzmann framework fails. These include: self-gravitating systemsPlastino and Plastino (1993); Aly (1993), turbulenceArimitsu and Arimitsu (2000a, b); Beck (2000); Beck et al. (2001), anomalous diffusionPlastino and Plastino (1995); Tsallis and Bukman (1996); Bologna et al. (2000); Malacarne et al. (2001), velocities of galaxiesLavagno et al. (1998), solar neutrinosLavagno and Quarati (2000), etc.

In spite of these successes, some drawbacks were identified in the early formalism, namely: (a) the density operator was not invariant under a uniform translation of energy spectrum; (b) the $`q`$-expected value of the identity operator was not the unity; (c) energy was not conservedTsallis et al. (1998). In a more recent work, Tsallis et al.Tsallis et al. (1998) circumvent these difficulties introducing the *normalized* constraint to the internal energy of the system.

Concerning applications to magnetic systems, the unormalized formalism was used by Portesi et al.Portesi et al. (1995) and Nobre et al.Nobre and Tsallis (1996), to describe the paramagnetic behavior of a system with N spins $`1/2`$. The authors found a non-measurable magnetic susceptibility, since it was exponentially dependent on the number N of particles. Mart nez et al.Mart nez et al. (2000) analyzed the same S=1/2 system in the framework of the normalized formalismTsallis et al. (1998), finding a similar result to that reported by the Portesi et al. and Nobre et al.

In the present work, we analyzed the paramagnetic behavior of N spins S, within the normalized formalism and showed that the system effective temperature T does not relate with the inverse of the Lagrange multiplier $`\beta `$, as it is normally assumed, but with the inverse of a re-scaled parameter $`\beta ^{}`$ ($`T=1/k\beta ^{}`$), already introduced by TsallisTsallis et al. (1998). In this approach, the unreal result found by PortesiPortesi et al. (1995), NobreNobre and Tsallis (1996) and Mart nezMart nez et al. (2000) was no longer found and the susceptibility became proportional to number of particles, as it is experimentally expected. In these early works, no attempt was made to correlate theoretical results to experimental magnetic systems of any kind.

We also suggest that the manganites are physical systems where the present concepts can be tested. Experimental data of Amaral et al.Amaral et al. (2001) and H bert et al.H bert et al. (2002), are in qualitative agreement with the results here reported. An analysis of non-extensivity in manganites in the ferromagnetic phase, also taking $`T=1/k\beta ^{}`$ as an effective temperature, was already published elsewhereReis et al. (2002).

## II The Model

We consider the Hamiltonian of a single spin S in a static and homogeneous magnetic field:

$$\widehat{}=\widehat{\stackrel{}{\mu }}\stackrel{}{B}=\widehat{\mu }_zB=g\mu __B\widehat{S}_zB$$
(1)
where g=2, for J=S. The $`q`$-generalized magnetic moment thermal average $`\widehat{\mu }_z_q`$ is:

$$\widehat{\mu }_z_q=g\mu _\text{B}\widehat{S}_z_q=g\mu _\text{B}SB_\text{S}^{\text{(q)}}$$
(2)
where $`\widehat{S}_z_q`$ is $`q`$-generalized thermal average spin operator and $`B_\text{S}^{\text{(q)}}`$ is the *generalized Brillouin function*. This first quantity can be determined in the Tsallis framework of normalized $`q`$-expectation valuesTsallis et al. (1998), from:

$$\widehat{S}_z_q=\frac{Tr\{\widehat{\rho }^q\widehat{S}_z\}}{Tr\{\widehat{\rho }^q\}}$$
(3)
where $`\widehat{\rho }`$ is the density operator, derived from the maximization of the entropy. Below, we analyze two different proposals for the entropy and show that re-scaling the Lagrange parameter, the same density operator emerges from both definitions.

### II.1 Tsallis Entropy

The entropy of the system is defined asTsallis (1988); Curado and Tsallis (1991); Tsallis et al. (1998):

$$𝒮_q=Tr\{\widehat{\rho }^q\widehat{𝒮}_q\}=\frac{k}{(q1)}[1Tr\{\widehat{\rho }^q\}]$$
(4)
where,

$$\widehat{𝒮}_q=k\mathrm{ln}_q\widehat{\rho }$$
(5)
and $`k`$ a positive constant. The generalized logarithm is defined asBorges (1998):

$$\mathrm{ln}_qf=(1q)^1(f^{1q}1)$$
(6)
The corresponding density operator $`\widehat{\rho }`$ can be determined from the maximization of the $`𝒮_q`$ functional, subjected to $`q`$-normalized constraintTsallis et al. (1998):

$$U_q=\frac{Tr\{\widehat{\rho }^q\widehat{}\}}{Tr\{\widehat{\rho }^q\}}$$
(7)
and the normalization of density operator $`Tr\{\widehat{\rho }\}=1`$, leading to:

$$\widehat{\rho }=\frac{1}{Z_q}\left[1(1q)\frac{\beta }{Tr\{\widehat{\rho }^q\}}(\widehat{}U_q)\right]^{\frac{1}{1q}}$$
(8)
where $`Z_q`$ is the generalized partition function:

$$Z_q=Tr\left\{\left[1(1q)\frac{\beta }{Tr\{\widehat{\rho }^q\}}(\widehat{}U_q)\right]^{\frac{1}{1q}}\right\}$$
(9)
and $`\beta `$ is the Lagrange parameter associated to internal energy constraintTsallis et al. (1998); Tsallis (2001). In this case, the entropy has a well defined concavity, for any value of $`q`$, being concave for $`q>0`$, and convex for $`q<0`$Tsallis (1988); Curado and Tsallis (1991); Tsallis et al. (1998); Tsallis (2001).

### II.2 Normalized Tsallis Entropy

If the entropy functional $`\widehat{𝒮}_q`$ is defined as proposed by Rajagopal and AbeRajagopal and Abe (1999):

$$𝒮_q=\frac{Tr\{\widehat{\rho }^q\widehat{𝒮}_q\}}{Tr\{\widehat{\rho }^q\}}=\frac{k}{(q1)}\left[\frac{1}{Tr\{\widehat{\rho }^q\}}1\right]$$
(10)
the density operator $`\widehat{\rho }`$ can be determined in a similar way as above, yielding:

$$\widehat{\rho }=\frac{1}{Z_q}\left[1(1q)\beta Tr\{\widehat{\rho }^q\}(\widehat{}U_q)\right]^{\frac{1}{1q}}$$
(11)
$$Z_q=Tr\left\{\left[1(1q)\beta Tr\{\widehat{\rho }^q\}(\widehat{}U_q)\right]^{\frac{1}{1q}}\right\}$$
(12)
In this case, the $`q`$ parameter is restricted to the interval $`0q1`$, preserving the entropy concavityRajagopal and Abe (1999).

The density operators emerging from both cenarios (Eq.8 and 11) can be re-written in the form:

$$\widehat{\rho }=\frac{1}{Z_q^{}}\left[1(1q)\beta ^{}\widehat{}\right]^{\frac{1}{1q}}$$
(13)
$$Z_q^{}=Tr\left\{\left[1(1q)\beta ^{}\widehat{}\right]^{\frac{1}{1q}}\right\}$$
(14)
where,

$$\beta ^{}=\frac{\beta }{Tr\{\widehat{\rho }^q\}+(1q)U_q\beta }$$
(15)
for the case A, and:

$$\beta ^{}=\frac{\beta }{\frac{1}{Tr\{\widehat{\rho }^q\}}+(1q)U_q\beta }$$
(16)
for case B. Here, we suggested that the effective temperature is:

$$T=\frac{1}{k\beta ^{}}$$
(17)
and the density operator becomes independent of the initial entropy functional.

The magnetic behavior of a S-spin system will be analyzed as a function of the parameter $`x^{}=g\mu _\text{B}SB\beta ^{}`$. In fact, in terms of $`x^{}`$, the generalized Brillouin function, is given by:

$$B_S^{(q)}=\frac{1}{S}\widehat{S}_z_q=\frac{1}{S}\frac{\underset{m_s=S}{\overset{+S}{}}m_s\left[1+(1q)x^{}\frac{m_s}{S}\right]^{\frac{q}{1q}}}{\underset{m_s=S}{\overset{+S}{}}\left[1+(1q)x^{}\frac{m_s}{S}\right]^{\frac{q}{1q}}}$$
(18)
It is to be remarked that cut-off procedureTsallis (2001); Lima and Penna (1999); Plastino and Plastino (1999); Tsallis (1999) implies that those states that do not satisfy the condition:

$$1+(1q)x^{}\frac{m_s}{S}0$$
(19)
must be excluded from the summation. In other words, these states are assigned with a zero probability amplitude, preserving the positive definition of the density operatorPortesi et al. (1995).

For a system constituted by N spins S particles the $`q`$-generalized spin operator thermal average (Eq.3), can be written as:

$$\widehat{S}_z_{q,N}=\frac{\underset{m_s=NS}{\overset{+NS}{}}Y(m_s)m_s\left[1+(1q)x^{}\frac{m_s}{S}\right]^{\frac{q}{1q}}}{\underset{m_s=NS}{\overset{+NS}{}}Y(m_s)\left[1+(1q)x^{}\frac{m_s}{S}\right]^{\frac{q}{1q}}}$$
(20)
where, $`Y(m_s)`$ is the multiplicity and :

$$\underset{m_s=NS}{\overset{NS}{}}Y(m_s)=(2S+1)^N$$
(21)
In the particular case of N $`1/2`$-spins particles, $`Y(m_s)`$ is simply given by:

$$Y(m_s)=\frac{N!}{\left(\frac{N}{2}m_s\right)!\left(\frac{N}{2}+m_s\right)!}$$
(22)
## III Results and Discussion

When $`q`$ is different from unity, general expressions for magnetic observables are difficult to calculate and interpret, even for simple systems. This hinders the comparison between theoretical predictions and experimental results. Numerical methods, on the other hand, provide a means to directly calculate observables, allowing a deeper understating of the theoretical results.

Figure 1 displays the generalized Brillouin function $`B_\text{S}^{\text{(q)}}`$ vs. $`x^{}`$ (Eq.18), for different values of $`q`$ and S=5/2. For $`q`$ up to 0.5, a series of kinks appear in the curve. One also notes that the saturation value $`B_\text{S}^{\text{(q)}}|_{\text{SAT}}`$ decreases with decreasing $`q`$. This is more clearly shown in figure 2(a) and (b), where $`B_\text{S}^{\text{(q)}}|_{\text{SAT}}`$ is plotted as a function of $`q`$, for half-integer and integer spin values, respectively. The behavior of a classical spin is also included and can be exactly calculated from Eq.18 with $`S\mathrm{}`$, giving:

$$B_\text{S}^{\text{(q)}}|_{\text{SAT}}=\frac{1}{2q}$$
(23)
The occupation probability (OP), as a function of $`x^{}`$, for each energy level of S=5/2 are displayed in figure 3(a)(b)(c), for several values of $`q`$. Figure 3(d) shown the same quantity for S=2 and q=0.1. We observe that the OP does not vanish for negative energy levels ($`m_s>0`$), even for very large values of $`x^{}`$, in contrast to what happens in the case q=1 (Maxwell-Boltzmann). From figure 3(c) we can see that, for q=0.1, the OP of the positive energy states ($`m_s<0`$) vanish sharply at the same $`x^{}`$ values as the kinks observed in figure 1. This occurs for q\<0.5. Therefore, we correlate the kinks observed on the magnetization curves to the lost of occupation of the most energetic states. This is consequence of Tsallis’ cut-off. In fact, from Eq.19, the $`x^{}`$ values where the kinks occurs can be derived, following:

$$x_{m_s}^{kinks}=\frac{S}{|m_s|(1q)}$$
(24)
It is to be remarked that for a half-integer spin, S+1/2 kinks occur, and S for an integer spin.

Figure 4(a)(b) display, respectively, the unormalized (Eq.4) and normalized (Eq.10) entropy $`𝒮_q`$, for S=5/2 and different $`q`$ values. Note that, for q\<0.5, the kinks discussed before are present. Besides, for any q\<1, the entropy does not vanish in the limit $`x^^10`$. In other words, even for high field and/or low temperature, the system has a finite entropy, which prevents a fully magnetized state. These features are valid for both entropy functionals (section II).

A general expression for the magnetic susceptiblity can be deduced from Eq.3 and Eq.18:

$$\chi _q=\underset{B0}{lim}\left[\frac{\widehat{\mu }_z_q}{B}\right]=\frac{C^{(q)}}{T}$$
(25)
where $`C^{(q)}=qC^{(1)}`$ is the generalized Curie constant. Figure 5 displays $`\chi _q^1`$ as a function of $`x^^1`$ for q=1.0 and 0.8, with S=5/2.

In order to compare the predictions of the proposed model to experimental data, we must investigate how the magnetic observables discussed before scales upon increasing the number of particles in the system. An extensive quantity $``$ scales as

$$_N=N_1$$
(26)
where N is the system number of particle.

Particularly useful for the purpose of comparison, is the magnetic susceptibility. Portesi et al.Portesi et al. (1995) and Nobre et al.Nobre and Tsallis (1996), considered the generalized magnetization on paramagnetic fase of a N $`1/2`$-spin system in the unormalized formalismCurado and Tsallis (1991); Tsallis et al. (1998); Portesi et al. (1995); Nobre and Tsallis (1996); Tsallis (2001). They found for $`\chi _q`$:

$$\chi _q=\frac{C_{\text{S=1/2}}^{\text{(q)}}}{T}2^{N(1q)}$$
(27)
where,

$$C_{\text{S=1/2}}^{\text{(q)}}=\frac{(g\mu _\text{B})^2}{4k}Nq$$
(28)
and $`T=1/k\beta `$ is the temperature of the system. Here, $`\beta `$ is the usual Lagrange parameter associated to the non-normalized constraint of internal energy.

Therefore, as N increases, the magnetic susceptibility becomes infinity for q\<1 and zero for q\>1. In other words, there is no linearity among $`\chi _q`$ and N. This phenomenon is called *dark magnetism*, in analogy to the cosmological concept of *dark matter* \[Portesi et al. (1995) and references there in\].

S. Mart nez et al.Mart nez et al. (2000), analyzed the paramagnetic behavior of the same N $`1/2`$-spin system, however, within the normalized formalismTsallis et al. (1998); Tsallis (2001), and found a similar result as Portesi and Nobre for $`\chi _q`$.

Our proposal is that the paramagnetic behavior of a non-extensive N $`1/2`$-spin system should be analyzed in the normalized formalism, using the density operator described in Eq.13, taking $`\beta ^{}`$, instead $`\beta `$, inversely proportional to the system temperature T. By doing so, the generalized paramagnetic susceptibility becomes:

$$\chi _q=\frac{C_{\text{S=1/2}}^{\text{(q)}}}{T}$$
(29)
which does not diverge as N increases, and is indeed proportional to N. Besides, this result is independent of the choice for the functional entropy.

## IV Possible Connection to Experimental Results

Amaral et al.Amaral et al. (2001) discussed the magnetic behavior of manganese oxides, namely, La0.67Ca0.33MnO3 and verified steps on the curve of M-1 vs. T in the paramagnetic phase, as shown in figure 6(a). These steps are analogue to those shown on figure 6(b), that represents also the inverse of magnetization as a function of $`x^1`$, for q=0.1, and encourage the idea of manganites as non-extensive objectsReis et al. (2002). The authors relate the change in the slope of the curve as an indication of cluster formation, which change the effective moment of Mn ions. These clusters could give rise to fractal structures, as discussed by DagottoDagotto et al. (2001), and therefore are in accordance to the ideas discussed here.

In addition, H bert et al.H bert et al. (2002) found magnetization curves in Pr0.5Ca0.5Mn0.95Ga0.05O3 that are also in qualitative agreement to those curves in figure 1.

## V Conclusion

In summary, we investigated the properties of a paramagnetic S-spin system under Tsallis generalized statistics, on the normalized formalism. For q\<0.5, a series of kinks appears in magnetization and entropy. This effect is a direct consequence of a peculiar occupation probability, as a result of Tsallis’ cut-off, where the positive energy states ($`m_s<0`$) vanish sharply. Additionally, the negative energy states ($`m_s>0`$) share non-zero occupation probability, preventing a fully magnetized state and the saturation magnetization decreases with decreasing $`q`$. We present evidences based on experimental results of Amaral et al.Amaral et al. (2001) and H bert et al.H bert et al. (2002) which add and support our previous publicationReis et al. (2002), where manganites are suggested to be magnetically non-extensive objects.

###### Acknowledgements.

One of us, M.S.R., acknowledges CAPES, Brazil, for financial support on Aveiro University. I.S.O. acknowledges FAPERJ, the Rio de Janeiro Funding Agency.

## VI Captions

Figure 1: Generalized Brillouin function for several $`q`$ values and S=5/2, as a function of $`x^{}`$.

Figure 2: Saturation magnetization $`B_\text{S}^{\text{(q)}}|_{\text{SAT}}`$ as a $`q`$ function, for (a)half-integer and (b)integer S values.

Figure 3: Probability of occupation of each energy level, for (a)(b)(c) S=5/2 and (d) S=2, as a function of $`x^{}`$, and various values of $`q`$.

Figure4: (a) Unormalized (Eq. 4) and (b) normalized (Eq. 10) entropy $`𝒮_q`$, for S=5/2 and different $`q`$ values, as a function of $`x^^1`$.

Figure5: Inverse of susceptibility for (a) q=1.0 and 0.8, for S=5/2, as a $`x^^1`$ function.

Figure6: (a) Temperature dependence of the inverse susceptibility for the manganite La2/3Ca1/3MnO3, above the Curie temperature ( 267 K), at low magnetic field. (b)Inverse of magnetization (Eq.18) as a function of $`x^1`$, for q=0.1