text:
# Taming chaos by impurities in two-dimensional oscillator arrays

## I Introduction

Coupled arrays of oscillators are studied extensively in many fields of science because of their prevalence in nature. They are used as models for coupled arrays of neurons , chemical reactions , coupled lasers  or Josephson junctions , charge-density-wave conductors , crystal dislocations in metals , and proton conductivity in hydrogen-bonded chains . Various models and coupling schemes have been proposed and analyzed previously . A particular class are arrays of coupled oscillators, which exhibit chaotic motion when uncoupled. This class includes the forced Frenkel-Kontorova model , which finds a straightforward physical realization in an array of diffusively coupled Josephson junctions , in which the applied current of each junction is modulated by a common frequency. The possibility to obtain synchronized motion in such systems has been investigated recently by Braiman et al. for the case of one- (1D) and two-dimensional (2D) chaotic arrays of forced damped nonlinear pendula  and coupled Josephson junctions . They observed the emergence of complex but frequency-locked spatio-temporal patterns, in which the chaotic behavior was completely suppressed, when a certain amount of disorder had been introduced by randomizing the lengths of the pendula. In Ref.  the same phenomenon has been investigated from a completely different point of view: It was shown for 1D arrays of coupled chaotic pendula, that introducing a single impurity at a particular site is sufficient to lead to complete synchronization.

In this paper we study 2D arrays of coupled chaotic pendula. We ask for the minimal coupling as well as the influence of concentration and arrangement of impurities, needed to observe spatio-temporal patterns. Although for geometrical reasons one might expect that a single impurity cannot play the role it plays in 1D arrays, we find that it is able to tame the chaotic behavior of an arbitrarily large 2D array, provided that the coupling is strong enough. A single impurity can produce synchronized spatio-temporal patterns, even though all oscillators and the impurity are chaotic when uncoupled. In order to observe spatio-temporal patterns for smaller couplings, the geometrical arrangement of the impurities is shown to play a crucial role. Specifically we show that an impurity configuration that divides the lattice into at least two disjoint parts is most appropriate for the creation of synchronized solutions. Such configurations will always produce patterns that locally can be identified as lines, provided the coupling is above a threshold value. The resulting spatio-temporal patterns are stable with respect to an increase of the system size, indicating that an impurity concentration that tends to zero may suffice to lead to synchronization.

## II Model

We will focus our analysis on the model examined in Refs. :

$`l_{n,m}^2\ddot{\theta }_{n,m}+\gamma \dot{\theta }_{n,m}=gl_{n,m}\mathrm{sin}\theta _{n,m}+\tau ^{}+\tau \mathrm{sin}\omega t`$ (1)
$`+k(\theta _{n+1,m}+\theta _{n1,m}4\theta _{n,m}+\theta _{n,m+1}+\theta _{n,m1}),`$ (2)
where $`n,m=1,2\mathrm{}N`$. Thus, there is a damped, driven pendulum with unity mass and length $`l_{nm}`$ on each site $`(n,m)`$ of the lattice, subject to an ac and a dc torque. The parameters used are the gravitational acceleration $`g=1`$, the dc torque $`\tau ^{}=0.7155`$, the ac torque $`\tau =0.4`$, the angular frequency $`\omega =0.25`$, and the damping $`\gamma =0.75`$. Neighbouring pendula are coupled via a discrete Laplacian, where $`k`$ denotes the coupling strength. We have chosen free boundary conditions, i.e. $`\theta _{0,m}=\theta _{1,m}`$, $`\theta _{N,m}=\theta _{N+1,m},\theta _{n,0}=\theta _{n,1}`$, $`\theta _{n,N}=\theta _{n,N+1}`$ and used a fourth order Runge-Kutta routine with a time step $`dt=0.01`$ to numerically integrate Eq. (1). We carefully checked that decreasing the time step to $`dt=0.001`$ did not alter our results.

A very convenient measure that allows a quick visualization of the average global spatio-temporal behavior of the lattice, is the average velocity

$$\sigma (jT)=\frac{1}{NM}\underset{n=1}{\overset{NM}{}}\dot{\theta }_n(jT)$$
(3)
at multiple times of the driving period $`T=1/\omega `$. Using this quantity  to obtain a bifurcation diagram not only can ascertain if chaotic or periodic behavior is obtained, but in addition helps to identify the maximum period of a pattern: Computing $`\sigma (t)`$ at each period $`t=jT`$ of the driving, will lead to a periodic sequence $`\sigma _1,\mathrm{},\sigma _p,\sigma _1,\mathrm{}`$, when transients have died out and a spatio-temporal pattern of periodicity $`p`$ (’P$`p`$ attractor’) has emerged. In practice we inspect the last $`20`$ values of $`\sigma (jT)`$ with $`j=1,\mathrm{}170`$, so that transients have died out. Thus, P20 attractors or attractors of larger periodicity are not recognized as such but rather appear as chaotic attractors.

Applying this strategy to an isolated pendulum, a bifurcation analysis with respect to the pendulum length $`l`$ was performed in Ref. . This approach revealed that each isolated pendulum is chaotic for values $`l=1\pm 0.002`$ and that three more chaotic windows exist: Two narrow ones at $`l0.84`$, and $`l0.52`$, and a broad one for $`l<0.35`$. Thus we know whether a chosen pendulum with length $`l`$ is chaotic or not. Furthermore, for the above chosen parameter values it is known, that pendula with $`l_{n,m}>1`$ and $`l_{n,m}<1`$ show libration and rotation, respectively, apart from the windows, where chaotic motion appears.

## III Results and Discussion

The simplest configuration of a 2D lattice described by Eq. (1), is that of a single impurity in a sea of identical chaotic pendula with length $`l_{n,m}=1`$. We fixed the lattice size to be $`50\times 50`$, with the the impurity located at site $`(25,25)`$ and made a bifurcation analysis with respect to the impurity length $`l_{\mathrm{imp}}`$ for various values $`k`$ of the coupling. In Fig. 1 the obtained bifurcation diagrams for $`k=1,2,3,5`$ are plotted versus $`l_{\mathrm{imp}}`$. For convenience we normalized $`\sigma `$ to take on values in the unit interval $`[0,1]`$. Similarly to 1D arrays we find, that one impurity is able to organize the 2D arrays. However, our bifurcation analysis shows that the coupling constant has to be larger than $`k_{\mathrm{cr}}3`$, whereas in the 1D case $`k0.1`$ was sufficient  to produce spatio-temporal patterns. We would like to point out the appearance of a P4 attractor for $`k=5,l_{\mathrm{imp}}<0.15`$ in Fig. 1d. Here an array of chaotic pendula is synchronized by a single impurity, which itself is chaotic, when isolated.

In order to illustrate the occurence of this P4 pattern in a better way, we show in Fig. 2 a typical gray scale plot of the velocities $`\dot{\theta }_{n,m}`$ as a function of the lattice coordinates $`(n,m)`$. The stroboscopic snapshots are taken with a time lag of one period of the driving $`T`$; darker shading indicates higher velocity. In this example we used $`64\times 64`$ oscillators having $`l_{n,m}=1`$ and located an impurity of length $`l_{imp}=0.1`$ in the middle of the lattice at $`(32,32)`$. The coupling constant is $`k=5`$. After some transient (not shown) the synchronization is maintained and a P4 attractor is clearly visible. We neglect edge effects except to remark that they are strictly confined to the last few pendula near the boundaries. Moreover, they become very regular if the calculation is carried on for longer times.

In many applications, however, one is interested in the weak coupling limit . In this limit, our analysis shows that one impurity is not able to create spatio-temporal organization of a chaotic lattice. Therefore, it is meaningful to ask whether spatio-temporal patterns can emerge at all and under which conditions they may be obtained. To this end we increase the number of impurities introduced in the lattice and investigate the importance of their geometrical arrangement as well as their concentration.

In Fig. 3a we report the P1 pattern emerging for a $`128\times 128`$ lattice with coupling constant $`k=0.5`$, when $`128`$ impurities were arranged along a line in the middle of the lattice, dividing it into two disjoint parts. The same behavior could be observed even for smaller couplings. We verified this for coupling constants as small as $`k_{\mathrm{cr}}=0.1`$. Moreover, we found that increasing the size of the lattice to $`N=256`$ (the limit of our computational capability), but maintaining the line-like geometry of impurities and the parameter-values, did not affect the formation of a pattern. In all cases we observed synchronization to a P1 pattern after an initial transient. Thus, increasing the size of the lattice will not affect the pattern formation although the percentage of the impurities will tend to zero in the limit of infinite systems. This indicates that the concentration of impurities is not of primary importance.

To test the influence of the special arrangement of impurities, we considered a $`128\times 128`$ lattice with $`k=0.5`$ and $`128`$ impurities of length $`l_{\mathrm{imp}}=0.7`$ located at random positions of the lattice. In all cases we have studied we obtained chaotic patterns. Even an increase of their concentration by more than a factor of 3 did not produce any spatio-temporal pattern.

It is thus natural to ask wether a line of impurities is the only geometry that can produce spatio-temporal patterns for relatively small values of the coupling constant. In Fig. 3b,c we report the P1 attractors of a lattice of $`128\times 128`$ coupled pendula with lengths $`l_{n,m}=1`$, $`l_{imp}=0.7`$ and coupling $`k=0.5`$, where the impurities are arranged like a cross (Fig. 3b) or as a ring (Fig. 3c). In both cases the observed pattern geometries consist locally of stripes as was also observed for the ’line-geometry’ in Fig. 3a.

What is the common geometrical feature of the above impurity configurations that allow them to control the chaotic lattice? From the above analysis we draw the conclusion, that it is sufficient for synchronization, that the impurites divide the lattice into at least two disjoint subregions. In that way the critical coupling needed to observe a spatio-temporal pattern is decreased by an order of magnitude. Reducing the coupling below $`k_{\mathrm{cr}}0.1`$, however, suppresses the formation of spatio-temporal patterns even for these lattices. This finding is consistent with earlier investigations on 1D arrays, which revealed a critical coupling $`k0.1`$, below which no pattern formation could be observed. Furthermore, the observed P1 pattern for the ’line geometry’ (Fig. 3a) is analogous to the 1D case, where the introduction of a single impurity caused the same topology of the array, i.e. a division into two disoint sets, and a similar pattern was observed . Thus, the results of the 1D case define the limiting $`k`$-value also for the 2D lattice.

We finally examined the effect of disorder and the possibility of obtaining self-organization or frequency locking. For a $`128\times 128`$ lattice we randomly varied the lengths of the pendula but restricted the range of the disorder such, that each individual pendulum is chaotic, i.e. $`l_{n,m}[0.998,1.002]`$. We found that the emerging pattern was always chaotic for a large number of different realizations of disorder and initial conditions. The lacking synchronization can be observed in a better way by inspecting $`\sigma (jT)`$, which does not show any periodicity (see Fig. 4a for an example). If the range of disorder is increased to include also regular moving pendula, i.e. $`l_{n,m}[0.8,1.2]`$, self organization is possible in agreement with findings of Ref. . This can be understood by the observation that synchronization already occurs, when choosing the disorder from an interval of pendulum lengths associated with regular motion, i.e. $`l_{n,m}[0.8,1.2]\backslash [0.998,1.002]`$. A particular example for this is shown in Figs. 4b,c, where the occurence of a P6 pattern can be observed. Thus, taking the lengths from the entire interval $`l_{n,m}[0.8,1.2]`$ on average yields only 1% chaotic pendula, whose motion will be overdominated by the synchronizing regular ones, explaining the spatio-temporal pattern observed in Ref. .

## IV Conclusion

In conclusion we have demonstrated that a lattice of chaotic pendula can be frequency locked into a spatio-temporal pattern by introducing impurities in the lattice. In the strong coupling limit, a single impurity can tame chaos. Decreasing the coupling constant requires more impurities in order to observe self-organization. In this case the geometry of the impurity configuration plays an important role. Our results suggest that if the impurity configuration divides the lattice into at least two disjoint parts then the coupling constant may be decreased without affecting the synchronization of the chaotic array. Moreover, the induced spatio-temporal patterns are then unaffected by the size of the lattice. Below the critical coupling $`k_{\mathrm{cr}}0.1`$ no synchronization is observed. The value of this critical coupling is dictated by the minimum coupling which leads to the formation of spatio-temporal patterns in the 1D case.

We are grateful to T. Gavrielides, V. Kovanis, A. Politi, and G Tsironis for helpful discussions.


format: text

text:
# Universal Behavior of Load Distribution in Scale-free Networks

## Abstract

We study a problem of data packet transport in scale-free networks whose degree distribution follows a power-law with the exponent $`\gamma `$. Load, or betweenness centrality of a vertex is the accumulated total number of data packets passing through that vertex when every pair of vertices send and receive a data packet along the shortest path connecting the pair. It is found that the load distribution follows a power-law with the exponent $`\delta 2.2(1)`$, insensitive to different values of $`\gamma `$ in the range, $`2<\gamma 3`$, and different mean degrees, which is valid for both undirected and directed cases. Thus, we conjecture that the load exponent is a universal quantity to characterize scale-free networks.

PACS numbers: 89.70.+c, 89.75.-k, 05.10.-a

Complex systems consist of many constituents such as individuals, substrates, and companies in social, biological, and economic systems, respectively, showing cooperative phenomena between constituents through diverse interactions and adaptations to the pattern they create. Interactions may be described in terms of graphs, consisting of vertices and edges, where vertices (edges) represent the constituents (their interactions). This approach was initiated by Erdös and Rényi (ER). In the ER model, the number of vertices is fixed, while edges connecting one vertex to another occur randomly with certain probability. However, the ER model is too random to describe real complex systems. Recently, Watts and Strogatz (WS)  introduced a small-world network, where a fraction of edges on a regular lattice is rewired with probability $`p_{\mathrm{𝑊𝑆}}`$ to other vertices. More recently, Barabási and Albert (BA)  introduced an evolving network where the number of vertices $`N`$ increases linearly with time rather than fixed, and a newly introduced vertex is connected to $`m`$ already existing vertices, following the so-called preferential attachment (PA) rule. When the number of edges $`k`$ incident upon a vertex is called the degree of the vertex, the PA rule means that the probability for the new vertex to connect to an already existing vertex is proportional to the degree $`k`$ of the selected vertex. Then the degree distribution $`P_D(k)`$ follows a power-law $`P_D(k)k^\gamma `$ with $`\gamma =3`$ for the BA model, while for the ER and WS models, it follows a Poisson distribution. Networks whose degree distribution follows a power-law, called scale-free (SF) networks, are ubiquitous in real-world networks such as the world-wide web, the Internet, the citation network and the author collaboration network of scientific papers , and the metabolic networks in biological organisms. On the other hand, there also exist random networks such as the actor network whose degree distribution follows a power-law but has a sharp cut-off in its tail. Thus, it has been proposed that the degree distribution can be used to classify a variety of diverse real-world networks. In SF networks, one may wonder if the exponent $`\gamma `$ is universal in analogy with the theory of critical phenomena; however, the exponent $`\gamma `$ turns out to be sensitive to the detail of network structure. Thus, a universal quantity for SF networks is yet to be found. From a theoretical viewpoint, it is important to find a universal quantity for SF networks, which is a purpose of this Letter.

A common feature between the WS and SF networks would be the small-world property that the mean separation between two vertices, averaged over all pairs of vertices (called the diameter hereafter), is shorter than that of a regular lattice. The small-world property in SF networks results from the presence of a few vertices with high degree. In particular, the hub, the vertex whose degree is the largest, plays a dominant role in reducing the diameter of the system. Diameters of many complex networks in real world are small, allowing objects transmitted through the network such as neural spikes on neural network, or data packets on Internet, to travel from one vertex to another quickly along the shortest path. The shortest paths are indeed of relevance to network transport properties. When a data packet is sent from one vertex to another through SF networks such as Internet, it is efficient to take a road along the shortest paths between the two. Then vertices with higher degrees should be heavily loaded and jammed by lots of data packets passing along the shortest paths. To prevent such Internet traffic congestions, and allow data packets to travel in a free-flow state, one has to enhance the capacity, the rate of data transmission, of each vertex to the extent that the capacity of each vertex is large enough to handle appropriately defined “load”.

In this Letter, we define and study such a quantity, which we simply call load, to characterize the transport dynamics in SF networks. In fact, this quantity turns out to be equivalent to “betweenness centrality” which was introduced in a social network to quantify how much power is centralized to people in social networks . While it has been noted that the betweenness centrality has a long tail, here we focus our attention on its probability distribution for various SF networks with different degree exponents. Thus knowing the distribution of such quantity enables us to not only estimate the capacity of each vertex needed for a free-free state, but also understand the power distribution in social networks, which is another purpose of this Letter.

To be specific, we suppose that a data packet is sent from a vertex $`i`$ to $`j`$, for every ordered pair of vertices $`(i,j)`$. For a given pair $`(i,j)`$, it is transmitted along the shortest path between them. If there exist more than one shortest paths, the data packet would encounter one or more branching points. In this case, we assume that the data packet is divided evenly by the number of branches at each branching point as it travels. Then we define the load $`\mathrm{}_k`$ at a vertex $`k`$ as the total amount of data packets passing through that vertex $`k`$ when all pairs of vertices send and receive one unit of data packet between them. Here, we do not take into account the time delay of data transfer at each vertex or edge, so that all data are delivered in a unit time, regardless of the distance between any two vertices.

We find numerically that the load distribution $`P_L(\mathrm{})`$ follows a power-law $`P_L(\mathrm{})\mathrm{}^\delta `$. Moreover, the exponent $`\delta 2.2`$ we obtained is insensitive to the detail of the SF network structure as long as the degree exponent is in the range, $`2<\gamma 3`$. The SF networks we used do not permit rewiring process, and the number of vertices are linearly proportional to that of edges. When $`\gamma >3`$, $`\delta `$ increases as $`\gamma `$ increases, however. The universal behavior is valid for directed networks as well, when $`2<\{\gamma _{\mathrm{in}},\gamma _{\mathrm{out}}\}3`$. Since the degree exponents in most of real-world SF networks satisfy $`2<\gamma 3`$, the universal behavior is interesting.

We construct a couple of classes of undirected SF networks both in the static and evolving ways. Each class of networks include a control parameter, according to which the degree exponent is determined. First, we deal with the static case. There are $`N`$ vertices in the system from the beginning, which are indexed by an integer $`i`$ $`(i=1,\mathrm{},N)`$. We assign the weight $`p_i=i^\alpha `$ to each vertex, where $`\alpha `$ is a control parameter in $`[0,1)`$. Next, we select two different vertices $`(i,j)`$ with probabilities equal to the normalized weights, $`p_i/_kp_k`$ and $`p_j/_kp_k`$, respectively, and add an edge between them unless one exists already. This process is repeated until $`mN`$ edges are made in the system. Then the mean degree is $`2m`$. Since edges are connected to a vertex with frequency proportional to the weight of that vertex, the degree at that vertex is given as

$$\frac{k_i}{_jk_j}\frac{(1\alpha )}{N^{1\alpha }i^\alpha },$$
(1)
where $`_jk_j=2mN`$. Then it follows that the degree distribution follows the power-law, $`P_D(k)k^\gamma `$, where $`\gamma `$ is given by

$$\gamma =(1+\alpha )/\alpha .$$
(2)
Thus, adjusting the parameter $`\alpha `$ in \[0,1), we can obtain various values of the exponent $`\gamma `$ in the range, $`2<\gamma <\mathrm{}`$.

Once a SF network is constructed, we select an ordered pair of vertices $`(i,j)`$ on the network, and identify the shortest path(s) between them and measure the load on each vertex along the shortest path using the modified version of the breath-first search algorithm introduced by Newman . It is found numerically that the load $`\mathrm{}_i`$ at vertex $`i`$ follows the formula,

$$\frac{\mathrm{}_i}{_j\mathrm{}_j}\frac{1}{N^{1\beta }i^\beta },$$
(3)
with $`\beta =0.80(5)`$. This value of $`\beta `$ is insensitive to different values of the exponent $`\gamma `$ in the range, $`2<\gamma 3`$ as shown in the inset of Fig. 1. The total load, $`_j\mathrm{}_j`$ scales as $`N^2\mathrm{log}N`$. This is because there are $`N^2`$ pairs of vertices in the system and the sum of the load contributed by each pair of vertices is equal to the distance between the two vertices, which is proportional to $`\mathrm{log}N`$. Therefore, the load $`\mathrm{}_i`$ at a vertex $`i`$ is given as

$$\mathrm{}_i(N\mathrm{log}N)(N/i)^\beta .$$
(4)
From Eq.(4), it follows that the load distribution scales as $`P_L(\mathrm{})\mathrm{}^\delta `$, with $`\delta =1+1/\beta 2.2(1)`$, independent of $`\gamma `$ in the range, $`2<\gamma 3`$. Direct measure of $`P_L(\mathrm{})`$ also gives $`\delta 2.2(1)`$ as shown in Fig. 1. We also check $`\delta `$ for different mean degrees $`m=2,4`$ and 6, but obtain the same value, $`\delta 2.2(1)`$. Thus, we conclude that the exponent $`\delta `$ is a generic quantity for this network. Note that Eqs.(1) and (4) combined gives a scaling relation between the load and the degree for this network as

$$\mathrm{}k^{(\gamma 1)/(\delta 1)}.$$
(5)
Thus, when and only when $`\gamma =\delta `$, the load at each vertex is directly proportional to its degree. Otherwise, it scales nonlinearly. On the other hand, for $`\gamma >3`$, the exponent $`\delta `$ depends on the exponent $`\gamma `$ in a way that it increases as $`\gamma `$ increases. Eventually, the load distribution decays exponentially for $`\gamma =\mathrm{}`$ as shown in Fig. 1. Thus, the transport properties of the SF networks with $`\gamma >3`$ are fundamentally different from those with $`2<\gamma 3`$. This is probably due to the fact that for $`\gamma >3`$, the second moment of $`P_D(\mathrm{})`$ exists, while for $`\gamma 3`$, it does not.

We examine the system-size dependent behavior of the load at the hub, $`\mathrm{}_h`$ for the static model. According to Eq.(4), $`\mathrm{}_h`$ behaves as $`\mathrm{}_hN^{1.8}\mathrm{log}N`$ in the range, $`2<\gamma 3`$, while for $`\gamma >3`$, $`\mathrm{}_h`$ increases with $`N`$ but at a much slower rate than that for $`2<\gamma 3`$ as shown in Fig. 2. That implies that the shortest pathways between two vertices become diversified, and they do not necessarily pass through the hub for $`\gamma >3`$. That may be related to the result that epidemic threshold is null in the range $`2<\gamma 3`$, while it is finite for $`\gamma >3`$ in SF networks, because there exist many other shortest paths not passing through the hub for $`\gamma >3`$, so that the infection of the hub does not always lead to the infection of the entire system. Thus, epidemic threshold is finite for $`\gamma >3`$ .

Next, we generate other SF networks in an evolving way, using the methods proposed by Kumar et al., and by Dorogovtsev et al.. In these cases, we also find the same results as in the case of static models.

We also consider the case of directed SF network. The directed SF networks are generated following the static rule. In this case, we assign two weights $`p_i=i^{\alpha _{\mathrm{out}}}`$ and $`q_i=i^{\alpha _{\mathrm{in}}}`$ $`(i=1,\mathrm{},N)`$ to each vertex for outgoing and incoming edges, respectively. Both control parameters $`\alpha _{\mathrm{out}}`$ and $`\alpha _{\mathrm{in}}`$ are in the interval $`[0,1)`$. Then two different vertices $`(i,j)`$ are selected with probabilities, $`p_i/_kp_k`$ and $`q_j/_kq_k`$, respectively, and an edge from the vertex $`i`$ to $`j`$ is created with an arrow, $`ij`$. The SF networks generated in this way show the power-law in both outgoing and incoming degree distributions with the exponents $`\gamma _{\mathrm{out}}`$ and $`\gamma _{\mathrm{in}}`$, respectively. They are given as $`\gamma _{\mathrm{out}}=(1+\alpha _{\mathrm{out}})/\alpha _{\mathrm{out}}`$ and $`\gamma _{\mathrm{in}}=(1+\alpha _{\mathrm{in}})/\alpha _{\mathrm{in}}`$. Thus, choosing various values of $`\alpha _{\mathrm{out}}`$ and $`\alpha _{\mathrm{in}}`$, we can determine different exponents $`\gamma _{\mathrm{out}}`$ and $`\gamma _{\mathrm{in}}`$. Following the same steps as for the undirected case, we obtain the load distribution on the directed SF networks. The load exponent $`\delta `$ obtained is $`2.3(1)`$ as shown in Fig. 3, consistent with the one for the undirected case, also being independent of $`\gamma _{\mathrm{out}}`$ and $`\gamma _{\mathrm{in}}`$ in $`2<\{\gamma _{\mathrm{out}},\gamma _{\mathrm{in}}\}3`$. Therefore, we conjecture that the load exponent is a universal value for both the undirected and directed cases.

To see if such universal value of $`\delta `$ appears in real world network, we analyzed the co-authorship network, where nodes represent scientists and they are connected if they wrote a paper together. The data are collected in the field of the neuro-science, published in the period 1991-1998. This network is appropriate to test the load, $`i.e.`$, the betweenness centrality distribution, because it does not include rewiring process as it evolves, and its degree exponent $`\gamma 2.2`$ lies in the range $`2<\gamma 3`$. As shown in Fig. 4, the load distribution follows a power-law with the exponent $`\delta 2.2`$, in good agreement with the value obtained in the previous models.

We also check the load distribution for the case when data travel with constant speed, so that the time delay of data transfer is proportional to the distance between two vertices. We find that the time delay effect does not change the load distribution and the conclusion of this work. The reason of this result is that when the time delay is accounted, load at each vertex is reduced roughly by a factor $`\mathrm{log}N`$, proportional to the diameter, which is negligible compared with the load without the time delay estimated to be $`N^{1.8}\mathrm{log}N`$ in Eq.(4). Due to this small-world property, the universal behavior remains unchanged under the time delay of data transmission.

Finally, we mention the load distribution of the small-world network of WS which is not scale-free. It is found that its load distribution does not obey a power-law, but shows a combined behavior of two Poisson-type decays resulting from short-ranged and long-ranged connections, respectively, as shown in Fig. 5. We also find the average load, $`\overline{\mathrm{}}(p_{\mathrm{𝑊𝑆}})(1/N)_i\mathrm{}_i(p_{\mathrm{𝑊𝑆}})`$, as a function of the rewiring probability $`p_{\mathrm{𝑊𝑆}}`$ decays rapidly with increasing $`p_{\mathrm{𝑊𝑆}}`$, behaving similar to the diameter in the WS model, as shown in the inset of Fig. 5.

In conclusion, we have considered a problem of data packet transport on scale-free networks generated according to preferential attachment rules and introduced a physical quantity, load $`\{\mathrm{}_i\}`$ at each vertex. We found that the load distribution follows a power-law, $`P_L(\mathrm{})\mathrm{}^\delta `$, with the exponent $`\delta 2.2(1)`$, which turns out to be insensitive to the degree exponent in the range $`(2,3]`$ when rewiring process is not included and the networks are of unaccerelated growth. Moreover, it is also the same for both directed and undirected cases within our numerical uncertainties. Therefore, we conjecture that the load exponent is a generic quantity to characterize scale-free networks. The universal behavior we found may have interesting implications to the interplay of SF network structure and dynamics. For $`\gamma >3`$, however, the load exponent $`\delta `$ increases as the degree exponent $`\gamma `$ increases, and eventually the load distribution decays exponentially as $`\gamma \mathrm{}`$. It would be interesting to examine the robustness of the universal behavior of the load distribution under some modifications of generating rules for SF networks such as rewiring process and acceleration growth, which, however, is beyond the scope of the current study.

We would like to thank M.E.J. Newman for introducing Refs. , and H. Jeong for providing the data for the co-authorship network. This work is supported by grants No.2000-2-11200-002-3 from the BRP program of the KOSEF.


text:
# Secular interactions between inclined planets and a gaseous disk

## 1 Introduction

The interaction between a young planetary system and its protoplanetary disk likely plays an important role in determining the orbital properties of the planets. Such interactions may be important for understanding the observed orbital properties of extra-solar planets (e.g. Marcy et al. 1999). Resonances within a gaseous disk likely play a key role in determining planetary eccentricities and inclinations. Much of the previous work, starting with Goldreich & Tremaine (1980), has emphasized the effects of mean-motion disk resonances, which involve frequencies that are comparable to the orbital frequencies of the planets.

Another potentially important class of resonances occurs where there is a matching of the precession frequency of a test particle (e.g. an asteroid) with the frequency of one of the precessional modes of the planetary system. The frequencies involved are much lower than in the mean-motion case. At such a resonance, the motion of a test particle can be strongly driven by the planets, resulting in a high orbital inclination (for a nodal resonance) or eccentricity (for an apsidal resonance). For example, in the solar system there is an important secular resonance that occurs near 2 AU due to driving involving Jupiter and Saturn. This resonance is believed to be associated with the inner truncation of the asteroid belt (Tisserand 1882).

It was recognized by Ward, Colombo, & Franklin (1976) that these resonances must have swept across portions of the early solar system owing to the effects of the gaseous disk, even if the planets do not migrate. The reason is that even a minimum-mass solar nebula can have an important influence on the relevant precession rates. As the nebula disperses, the precession rates vary, along with the resonance locations (Ward 1981). Regions through which the resonances sweep may be driven into significantly eccentric or inclined orbits, such as are observed among the asteroids today (e.g. Nagasawa, Tanaka, & Ida 2000). The resonances can also have an important bearing on the conditions for planet formation by stirring the planetesimal disk.

However, collective effects can be important within the disk and may even result in the development of density waves or bending waves. Consequently, the response is less localized than would be indicated by test particles, as was recognized by Ward & Hahn (1998) and Tremaine (1998). These investigations concentrated on the collective effects of self-gravity in particle disks.

Secular interactions involving the gas disk of a young planetary system are potentially of importance, since much more mass is contained in the gas disk than in the planetesimal disk. The disk is a fluid body and is capable of wave motions through (at least) compressive, inertial, buoyancy, and self-gravitational forces. It is also capable of dissipating energy through a turbulent effective viscosity or through shocks. It is therefore important to understand the dynamical response of the disk to secular forcing and the implications for the evolution of the planetary orbits.

The purpose of this paper is to investigate the outcome of secular interactions in mutually inclined planet-disk systems. The relevant waves in the disk are long-wavelength bending waves with azimuthal wavenumber $`m=1`$. The equations governing such waves in a protostellar disk are fairly well established, while those for long-wavelength eccentric density waves have received less attention (although see Ogilvie 2001). The importance of understanding the global $`m=1`$ response of the disk to secular perturbations has been emphasized by Tremaine (1998).

The general outline of the paper is as follows. Sections 2–5 describe the physical model for nodal secular interactions and formulate the normal mode analysis. Sections 6 and 7 explore the limit of sufficiently slow modes, when the disk responds nearly rigidly. Section 8 describes how the effects of mean-motion resonances are included in the normal mode analysis. Sections 9 and 10 describe numerical results for Jupiter interacting with the solar nebula. Section 11 provides a simple scaling analysis based on a nearly rigid tilt model. Section 12 discusses the relationship of the current approach with that of earlier work involving WKB theory. Section 13 describes the analysis of a close-orbiting planet in the central hole of a disk. Section 14 discusses Jupiter and Saturn interacting with the solar nebula at the $`\nu _{16}`$ secular resonance.

## 2 Modeling secular interactions

Planets experience secular gravitational interactions on time-scales much longer than their orbital periods. On such time-scales the planets may be considered as continuous rings representing their average mass density. These rings are in general elliptical and mutually inclined. Gravitational interactions between the rings lead to apsidal and nodal precession of the orbits.

In the case of small eccentricities and inclinations, the secular evolution of a system of $`n`$ planets can be described in terms of normal modes of a linear dynamical system, $`n`$ modes for eccentricity and $`n`$ for inclination (e.g. Murray & Dermott 1999). The eigenvector of a mode describes the relative distribution of eccentricity or inclination among the planetary orbits, while the eigenfrequency is the rate at which the pattern precesses.

In typical situations involving mean-motion resonances (Lindblad or vertical resonances), if a particle resonance is located at a certain radius in a fluid disk, the disk responds by launching a wave that carries energy and angular momentum away from the resonance. Such waves are usually assumed implicitly to be damped through some dissipative process before reaching an edge of the disk from which they might otherwise reflect. The gravitational torque between the disk and the perturber can then be calculated from a standard formula evaluated at the location of the resonance (Goldreich & Tremaine 1979).

However, in the case of a secular resonance, the frequency may be so small that the wavelength would be comparable to (or may even exceed) the size of the disk. In that case the global response of the disk must be computed, including explicit dissipation and the correct boundary conditions. The exact location of the resonance then ceases to have great importance and the new possibility of a global secular resonance arises.

We focus on the case of inclination, so that the relevant waves in the disk are long-wavelength bending waves with azimuthal wavenumber $`m=1`$. As noted above, the equations governing such waves in a protostellar disk are fairly well established.

## 3 Basic equations

Let $`(r,\varphi ,z)`$ be cylindrical polar coordinates with origin at the central object, of mass $`M_{}`$. We consider a protostellar disk of semi-thickness $`H(r)`$, for which the turbulent viscosity parameter $`\alpha `$ satisfies the condition $`\alpha H/r`$. The linearized equations for bending waves in such a disk have been derived in several papers (Papaloizou & Lin 1995; Masset & Tagger 1996; Demianski & Ivanov 1997). We present them in the form (Lubow & Ogilvie 2000)

$$\mathrm{\Sigma }r^2\mathrm{\Omega }\frac{W}{t}=\frac{1}{r}\frac{𝒢}{r}+T,$$
(1)
$$\frac{𝒢}{t}+\left(\frac{\kappa ^2\mathrm{\Omega }^2}{\mathrm{\Omega }^2}\right)\frac{i\mathrm{\Omega }}{2}𝒢+\alpha \mathrm{\Omega }𝒢=\frac{r^3\mathrm{\Omega }^3}{4}\frac{W}{r}.$$
(2)
Here $`W(r,t)=l_x+il_y`$ is the complex tilt variable, which describes the warped shape of the disk. Essentially, the disk can be thought of as a continuum of concentric circular rings with radii $`r`$ and unit normal vectors $`𝒍(r,t)`$. Also $`𝒢(r,t)=G_x+iG_y`$ is the complex internal torque variable. In a warped disk there is a horizontal internal torque $`2\pi 𝑮(r,t)`$ that is mediated by horizontal motions that are proportional to the distance $`z`$ from the midplane. These stresses are responsible for the propagation of bending waves. Similarly $`T(r,t)=T_x+iT_y`$ is the complex horizontal external torque density acting on the disk. Finally $`\mathrm{\Sigma }(r)`$ is the surface density, $`\mathrm{\Omega }(r)`$ the orbital angular velocity, $`\kappa (r)`$ the epicyclic frequency and $`(r)`$ the second vertical moment of the density. These are defined by

$$\mathrm{\Sigma }=\rho 𝑑z,=\rho z^2𝑑z,$$
(3)
$$r\mathrm{\Omega }^2=\frac{\mathrm{\Phi }}{r}|_{z=0},\kappa ^2=4\mathrm{\Omega }^2+2r\mathrm{\Omega }\frac{d\mathrm{\Omega }}{dr},$$
(4)
where $`\mathrm{\Phi }(r,z)`$ is the (axisymmetric component of the) gravitational potential experienced by the disk.

Consider a thin, uniform circular ring of mass $`m_iM_{}`$ and radius $`r_i`$, representing either an annulus of the disk or the time-average of a planetary orbit. Its contribution to the potential in the plane $`z=0`$ is

$$\mathrm{\Phi }_i=\frac{Gm_i}{2\pi }_0^{2\pi }(r^2+r_i^22rr_i\mathrm{cos}\varphi )^{1/2}𝑑\varphi .$$
(5)
Therefore the angular velocity and epicyclic frequency experienced by the ring are obtained from summations over all other such rings, in the forms

$$\mathrm{\Omega }_i^2=\frac{GM_{}}{r_i^3}+\underset{ji}{}\frac{2Gm_j}{r_ir_j}\left[K_0(r_i,r_j)\frac{r_j}{r_i}K_1(r_i,r_j)\right],$$
(6)
$$\kappa _i^2=\frac{GM_{}}{r_i^3}+\underset{ji}{}\frac{2Gm_j}{r_ir_j}\left[K_0(r_i,r_j)\frac{2r_j}{r_i}K_1(r_i,r_j)\right],$$
(7)
where $`K_0`$ and $`K_1`$ are symmetric kernels with dimensions of inverse length, given by

$$K_m(r_i,r_j)=\frac{r_ir_j}{4\pi }_0^{2\pi }(r_i^2+r_j^22r_ir_j\mathrm{cos}\varphi )^{3/2}\mathrm{cos}^m\varphi d\varphi .$$
(8)
Hereafter we write $`K_1`$ simply as $`K`$.

Now consider the tilt interaction between two such rings of masses $`m_i`$ and $`m_j`$, radii $`r_i`$ and $`r_j`$, and tilt vectors $`𝒍_i`$ and $`𝒍_j`$. For small relative inclinations, the gravitational torque exerted by ring $`j`$ on ring $`i`$ is

$$𝑻_{ji}=Gm_im_jK(r_i,r_j)𝒍_i\times 𝒍_j.$$
(9)
This expression can be obtained by averaging the torque exerted by a point mass on a circular ring (Lubow & Ogilvie 2000) over the orbital motion of the point mass. Similar expressions can be found in studies of galactic warps (e.g. Sparke & Casertano 1988). In the complex notation, for small inclinations from the $`xy`$-plane, the total external torque acting on ring $`i`$ is

$$T_i=\underset{ji}{}Gm_im_jK(r_i,r_j)i(W_jW_i).$$
(10)
The kernels can also be written as

$$K_m(r_i,r_j)=\frac{r_<}{4r_>^2}b_{3/2}^{(m)}\left(\frac{r_<}{r_>}\right),$$
(11)
where $`r_>=\mathrm{max}(r_i,r_j)`$, $`r_<=\mathrm{min}(r_i,r_j)`$ and $`b_\gamma ^{(m)}`$ is the Laplace coefficient. In practice we evaluate these in terms of elliptic integrals using Carlson’s algorithms (Press et al. 1992).

These kernels diverge as $`|r_ir_j|0`$. In reality the gravitational interaction remains bounded because of the non-zero vertical thickness of the rings, which we have neglected. To take account of this in an approximate way, we soften the kernels by replacing

$$K_m(r_i,r_j)K_m(\frac{r_i+r_jh}{2},\frac{r_i+r_j+h}{2}),$$
(12)
whenever $`|r_ir_j|<h`$. Here the smoothing length $`h`$ is to be understood as an approximate measure of the thickness of the rings.

## 4 Coupled systems of planets and disks

We consider a general system consisting of multiple planets embedded in a gaseous disk. Provided the planets create gaps in the disk, the disk then becomes partitioned into a set of disks. The disks and planets interact through gravity.

We consider planets $`i=1,2,\mathrm{},n_\mathrm{p}`$ of masses $`m_{\mathrm{p}i}`$ in circular orbits of radii $`r_{\mathrm{p}i}`$ and angular velocities $`\mathrm{\Omega }_{\mathrm{p}i}`$. We also consider disks $`k=1,2,\mathrm{},n_\mathrm{d}`$ of inner radii $`a_k`$ and outer radii $`b_k`$. The inclinations of the planetary orbits are described by the tilt variables $`W_{\mathrm{p}i}(t)`$. The tilt within a disk is described by $`W_{\mathrm{d}k}(r,t)`$. We assume throughout that the inclinations are small, so that linear theory applies.

If only secular nodal interactions are considered, the planets can be treated as inclined circular rings interacting with each other and with the disks through the torques described in Section 3. We return later (Section 8) to the effect of mean-motion resonances on tilt evolution. However, we assume throughout that any evolution of the surface density of the disk, or of the semi-major axes of the planetary orbits, may be neglected.

The secular equations for the planets are then of the form

$`m_{\mathrm{p}i}r_{\mathrm{p}i}^2\mathrm{\Omega }_{\mathrm{p}i}{\displaystyle \frac{dW_{\mathrm{p}i}}{dt}}`$ $`=`$ $`{\displaystyle \underset{ji}{}}Gm_{\mathrm{p}i}m_{\mathrm{p}j}K(r_{\mathrm{p}i},r_{\mathrm{p}j})i(W_{\mathrm{p}j}W_{\mathrm{p}i})`$ (13)
$`+{\displaystyle \underset{k}{}}{\displaystyle _{a_k}^{b_k}}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})i(W_{\mathrm{d}k}W_{\mathrm{p}i})\mathrm{\hspace{0.17em}2}\pi r𝑑r,`$
while the equations for the disks are

$`\mathrm{\Sigma }_kr^2\mathrm{\Omega }{\displaystyle \frac{W_{\mathrm{d}k}}{t}}`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{𝒢_k}{r}}+{\displaystyle \underset{i}{}}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})i(W_{\mathrm{p}i}W_{\mathrm{d}k})`$ (14)
$`+{\displaystyle \underset{\mathrm{}}{}}{\displaystyle _a_{\mathrm{}}^b_{\mathrm{}}}G\mathrm{\Sigma }_k\mathrm{\Sigma }_{\mathrm{}}^{}K(r,r^{})i(W_\mathrm{d}\mathrm{}^{}W_{\mathrm{d}k})\mathrm{\hspace{0.17em}2}\pi r^{}𝑑r^{},`$
$$\frac{𝒢_k}{t}+\left(\frac{\kappa ^2\mathrm{\Omega }^2}{\mathrm{\Omega }^2}\right)\frac{i\mathrm{\Omega }}{2}𝒢_k+\alpha \mathrm{\Omega }𝒢_k=\frac{_kr^3\mathrm{\Omega }^3}{4}\frac{W_{\mathrm{d}k}}{r},$$
(15)
where $`\mathrm{\Sigma }_{\mathrm{}}^{}=\mathrm{\Sigma }_{\mathrm{}}(r^{})`$ and $`W_\mathrm{d}\mathrm{}^{}=W_\mathrm{d}\mathrm{}(r^{},t)`$. The final term in equation (14) represents the self-gravitation of each disk and the gravitational interactions between disks.

The (complex) total horizontal angular momentum of the system,

$$\underset{i}{}m_{\mathrm{p}i}r_{\mathrm{p}i}^2\mathrm{\Omega }_{\mathrm{p}i}W_{\mathrm{p}i}+\underset{k}{}_{a_k}^{b_k}\mathrm{\Sigma }_kr^2\mathrm{\Omega }W_{\mathrm{d}k}\mathrm{\hspace{0.17em}2}\pi r𝑑r,$$
(16)
is exactly conserved if, as we assume, the boundary conditions

$$𝒢_k(a_k,t)=𝒢_k(b_k,t)=0$$
(17)
hold at all times $`t`$. There is also a conservation law for the vertical angular momentum, which we derive in Appendix A.

## 5 Normal modes

We have a set of coupled integro-differential equations that are linear and homogeneous. Solutions may be sought in the form of normal modes,

$$W_{\mathrm{p}i}(t)=\stackrel{~}{W}_{\mathrm{p}i}e^{i\omega t},W_{\mathrm{d}k}(r,t)=\stackrel{~}{W}_{\mathrm{d}k}(r)e^{i\omega t},𝒢_k(r,t)=\stackrel{~}{𝒢}_k(r)e^{i\omega t},$$
(18)
where $`\omega `$ is a complex frequency eigenvalue. The real part of $`\omega `$ is the precession rate of the mode, while the imaginary part is the decay rate. Substituting this and omitting the tildes, we obtain the integral and integro-differential equations

$`i\omega m_{\mathrm{p}i}r_{\mathrm{p}i}^2\mathrm{\Omega }_{\mathrm{p}i}W_{\mathrm{p}i}`$ $`=`$ $`{\displaystyle \underset{ji}{}}Gm_{\mathrm{p}i}m_{\mathrm{p}j}K(r_{\mathrm{p}i},r_{\mathrm{p}j})i(W_{\mathrm{p}j}W_{\mathrm{p}i})`$ (19)
$`+{\displaystyle \underset{k}{}}{\displaystyle _{a_k}^{b_k}}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})i(W_{\mathrm{d}k}W_{\mathrm{p}i})\mathrm{\hspace{0.17em}2}\pi r𝑑r,`$
$`i\omega \mathrm{\Sigma }_kr^2\mathrm{\Omega }W_{\mathrm{d}k}`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{d𝒢_k}{dr}}+{\displaystyle \underset{i}{}}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})i(W_{\mathrm{p}i}W_{\mathrm{d}k})`$ (20)
$`+{\displaystyle \underset{\mathrm{}}{}}{\displaystyle _a_{\mathrm{}}^b_{\mathrm{}}}G\mathrm{\Sigma }_k\mathrm{\Sigma }_{\mathrm{}}^{}K(r,r^{})i(W_\mathrm{d}\mathrm{}^{}W_{\mathrm{d}k})\mathrm{\hspace{0.17em}2}\pi r^{}𝑑r^{},`$
$$i\omega 𝒢_k+\left(\frac{\kappa ^2\mathrm{\Omega }^2}{\mathrm{\Omega }^2}\right)\frac{i\mathrm{\Omega }}{2}𝒢_k+\alpha \mathrm{\Omega }𝒢_k=\frac{_kr^3\mathrm{\Omega }^3}{4}\frac{dW_{\mathrm{d}k}}{dr},$$
(21)
subject to the boundary conditions

$$𝒢_k(a_k)=𝒢_k(b_k)=0.$$
(22)
Frequencies $`\mathrm{\Omega }`$ and $`\kappa `$ are given by equations (6) and (7). All orbits are assumed to be prograde. There is always a trivial rigid-tilt mode with $`\omega =0`$, $`𝒢_k=0`$, and $`W_{\mathrm{p}i}=W_{\mathrm{d}k}=\mathrm{constant}`$.

## 6 Rigid response

Secular modes usually have very long periods. If a disk is easily able to maintain radial communication (through pressure, viscosity or self-gravitation) on this long time-scale, it will participate in such a mode almost as a rigid body, and $`W`$ will not vary much across the disk.

Suppose that each disk indeed behaves rigidly ($`W_{\mathrm{d}k}=\mathrm{constant}`$). Then $`𝒢_k`$ can be eliminated by multiplying equation (20) by $`2\pi r`$, integrating from $`a_k`$ to $`b_k`$, and using the boundary conditions. We then find

$$\omega J_{\mathrm{p}i}W_{\mathrm{p}i}=\underset{ji}{}C_{ij}^{\mathrm{pp}}(W_{\mathrm{p}j}W_{\mathrm{p}i})+\underset{k}{}C_{ik}^{\mathrm{pd}}(W_{\mathrm{d}k}W_{\mathrm{p}i}),$$
(23)
$$\omega J_{\mathrm{d}k}W_{\mathrm{d}k}=\underset{i}{}C_{ik}^{\mathrm{pd}}(W_{\mathrm{p}i}W_{\mathrm{d}k})+\underset{\mathrm{}k}{}C_k\mathrm{}^{\mathrm{dd}}(W_\mathrm{d}\mathrm{}W_{\mathrm{d}k}),$$
(24)
where

$$J_{\mathrm{p}i}=m_{\mathrm{p}i}r_{\mathrm{p}i}^2\mathrm{\Omega }_{\mathrm{p}i}$$
(25)
is the angular momentum of planet $`i`$, and

$$J_{\mathrm{d}k}=_{a_k}^{b_k}\mathrm{\Sigma }r^2\mathrm{\Omega }\mathrm{\hspace{0.17em}2}\pi r𝑑r$$
(26)
is the total angular momentum of disk $`k`$. The coupling coefficients are defined by

$$C_{ij}^{\mathrm{pp}}=Gm_{\mathrm{p}i}m_{\mathrm{p}j}K(r_{\mathrm{p}i},r_{\mathrm{p}j}),$$
(27)
$$C_{ik}^{\mathrm{pd}}=_{a_k}^{b_k}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})\mathrm{\hspace{0.17em}2}\pi r𝑑r,$$
(28)
$$C_k\mathrm{}^{\mathrm{dd}}=_a_{\mathrm{}}^b_{\mathrm{}}_{a_k}^{b_k}G\mathrm{\Sigma }_k\mathrm{\Sigma }_{\mathrm{}}^{}K(r,r^{})\mathrm{\hspace{0.17em}2}\pi r𝑑r\mathrm{\hspace{0.17em}2}\pi r^{}𝑑r^{},$$
(29)
Evidently the disks can be treated formally as additional planets if they behave as rigid bodies. By introducing Greek indices that run from $`1`$ to $`n_\mathrm{p}+n_\mathrm{d}`$, we obtain the algebraic eigenvalue problem

$$\omega J_\beta W_\beta =\underset{\gamma \beta }{}C_{\beta \gamma }(W_\gamma W_\beta ).$$
(30)
This is formally identical to the multiple-planet problem described by, e.g., Murray & Dermott (1999). Since $`C`$ is a symmetric matrix, we have

$$\omega \underset{\beta }{}J_\beta |W_\beta |^2=\frac{1}{2}\underset{\beta ,\gamma }{}C_{\beta \gamma }|W_\beta W_\gamma |^2.$$
(31)
Since, further, the coefficients $`C_{\beta \gamma }`$ are real and positive, the eigenvalues are all real and non-positive. Apart from the trivial rigid-tilt mode, all modes precess retrogradely without growth or decay.

Suppose the system contains only two components, e.g. one planet interacting with a connected disk. The rigid solution satisfies

$$\left[\begin{array}{cc}\omega J_1+C& C\\ C& \omega J_2+C\end{array}\right]\left[\begin{array}{c}W_1\\ W_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right],$$
(32)
where $`C=C_{12}`$ is the coupling coefficient between the planet and the disk. The eigenvalues and eigenvectors are

$$\omega =0,\left[\begin{array}{c}W_1\\ W_2\end{array}\right]\left[\begin{array}{c}1\\ 1\end{array}\right],$$
(33)
a trivial solution corresponding to a rigid tilt, and

$$\omega =\frac{C(J_1+J_2)}{J_1J_2},\left[\begin{array}{c}W_1\\ W_2\end{array}\right]\left[\begin{array}{c}J_2\\ J_1\end{array}\right],$$
(34)
corresponding to a retrogradely precessing mode in which the tilt of each component is inversely proportional to its angular momentum (in order to conserve the total angular momentum).

In a three-component system, e.g. a planet interacting with interior and exterior disks, the rigid solutions satisfy

$$\left[\begin{array}{ccc}\omega J_1+C_{12}+C_{13}& C_{12}& C_{13}\\ C_{12}& \omega J_2+C_{12}+C_{23}& C_{23}\\ C_{13}& C_{23}& \omega J_3+C_{13}+C_{23}\end{array}\right]\left[\begin{array}{c}W_1\\ W_2\\ W_3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right],$$
(35)
where subscript 1 refers to the planet, and subscripts 2 and 3 to the inner and outer disks.

In the general case, the eigenvalues and eigenfunctions are algebraically complicated. However, consider the case in which the planet is much less massive than either disk. The first non-trivial mode is the equivalent of equation (34), but involving the two dominant components:

$$\omega \frac{C_{23}(J_2+J_3)}{J_2J_3},\left[\begin{array}{c}W_2\\ W_3\end{array}\right]\left[\begin{array}{c}J_3\\ J_2\end{array}\right].$$
(36)
The planet’s tilt is then driven according to

$$W_1=\frac{C_{12}W_2+C_{13}W_3}{\omega J_1+C_{12}+C_{13}}.$$
(37)
Secular resonance occurs here if the frequency is such that the denominator vanishes.

Secular resonance of this type could also occur, for example, in a system of two massive planets and one low-mass disk. In this case the disk would be responding globally as a rigid body. This is quite different from the situation in which a secular particle resonance is located somewhere within a disk.

The second non-trivial mode has

$$\omega \frac{(C_{12}+C_{13})}{J_1},|W_2|,|W_3||W_1|.$$
(38)
Here only the planet’s orbit is significantly tilted, and precesses at a rate determined by the two disks.

## 7 Nearly rigid response

The internal torque $`𝒢_k`$ required for a disk to respond rigidly can be determined from equation (20). In general, equation (21) will then indicate that a slight warping is in fact required. The rigid solution can be understood as the leading term in an expansion of the solution in powers of a small parameter. If $`ϵ=H/r`$ is the angular semi-thickness of the disk, and $`c_\mathrm{s}H\mathrm{\Omega }`$ the isothermal sound speed, the small parameter here is a characteristic value of $`(\omega /ϵ\mathrm{\Omega })^2`$, or $`(\omega r/c_\mathrm{s})^2`$. This parameter is small when the disk can communicate effectively, i.e. when the time taken for a bending wave to propagate across the disk is less than the wave period.

Now suppose the disk behaves nearly rigidly, so that we may pose an expansion in powers of this parameter, of the form

$`\omega `$ $`=`$ $`\omega ^{(0)}+\omega ^{(1)}+\mathrm{},`$ (39)
$`W_{\mathrm{p}i}`$ $`=`$ $`W_{\mathrm{p}i}^{(0)}+W_{\mathrm{p}i}^{(1)}+\mathrm{},`$ (40)
$`W_{\mathrm{d}k}`$ $`=`$ $`W_{\mathrm{d}k}^{(0)}+W_{\mathrm{d}k}^{(1)}(r)+\mathrm{},`$ (41)
$`𝒢_k`$ $`=`$ $`𝒢_k^{(0)}(r)+𝒢_k^{(1)}(r)+\mathrm{}.`$ (42)
At leading order we obtain the rigid solution described in Section 6, denoted by the superscript $`(0)`$. At the next order we find

$`i\omega ^{(1)}J_{\mathrm{p}i}W_{\mathrm{p}i}^{(0)}+i\omega ^{(0)}J_{\mathrm{p}i}W_{\mathrm{p}i}^{(1)}={\displaystyle \underset{ji}{}}Gm_{\mathrm{p}i}m_{\mathrm{p}j}K(r_{\mathrm{p}i},r_{\mathrm{p}j})i(W_{\mathrm{p}j}^{(1)}W_{\mathrm{p}i}^{(1)})`$
$`+{\displaystyle \underset{k}{}}{\displaystyle _{a_k}^{b_k}}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})i(W_{\mathrm{d}k}^{(1)}W_{\mathrm{p}i}^{(1)})\mathrm{\hspace{0.17em}2}\pi r𝑑r,`$ (43)
$`i\omega ^{(1)}\mathrm{\Sigma }_kr^2\mathrm{\Omega }W_{\mathrm{d}k}^{(0)}+i\omega ^{(0)}\mathrm{\Sigma }_kr^2\mathrm{\Omega }W_{\mathrm{d}k}^{(1)}={\displaystyle \frac{1}{r}}{\displaystyle \frac{d𝒢_k^{(1)}}{dr}}+{\displaystyle \underset{i}{}}Gm_{\mathrm{p}i}\mathrm{\Sigma }_kK(r,r_{\mathrm{p}i})i(W_{\mathrm{p}i}^{(1)}W_{\mathrm{d}k}^{(1)})`$
$`+{\displaystyle \underset{\mathrm{}}{}}{\displaystyle _a_{\mathrm{}}^b_{\mathrm{}}}G\mathrm{\Sigma }_k\mathrm{\Sigma }_{\mathrm{}}^{}K(r,r^{})i(W_\mathrm{d}\mathrm{}^{(1)}W_{\mathrm{d}k}^{(1)})\mathrm{\hspace{0.17em}2}\pi r^{}𝑑r^{},`$ (44)
$$i\omega ^{(0)}𝒢_k^{(0)}+\left(\frac{\kappa ^2\mathrm{\Omega }^2}{\mathrm{\Omega }^2}\right)\frac{i\mathrm{\Omega }}{2}𝒢_k^{(0)}+\alpha \mathrm{\Omega }𝒢_k^{(0)}=\frac{_kr^3\mathrm{\Omega }^3}{4}\frac{dW_{\mathrm{d}k}^{(1)}}{dr},$$
(45)
subject to the boundary conditions

$$𝒢_k^{(1)}(a_k)=𝒢_k^{(1)}(b_k)=0.$$
(46)
After some manipulations we may eliminate $`𝒢_k^{(1)}`$ by integration and obtain

$$\omega ^{(1)}\underset{\beta }{}J_\beta |W_\beta ^{(0)}|^2=\underset{k}{}_{a_k}^{b_k}\left(\frac{4}{_kr^4\mathrm{\Omega }^3}\right)\left[\omega ^{(0)}\left(\frac{\kappa ^2\mathrm{\Omega }^2}{\mathrm{\Omega }^2}\right)\frac{\mathrm{\Omega }}{2}+i\alpha \mathrm{\Omega }\right]|𝒢_k^{(0)}|^2\mathrm{\hspace{0.17em}2}\pi r𝑑r.$$
(47)
Since $`𝒢_k^{(0)}`$ does not depend on $`\alpha `$ but is proportional to $`\mathrm{\Sigma }`$, this shows that nearly rigid modes are damped at a rate $`\mathrm{Im}(\omega ^{(1)})`$ proportional to $`\alpha /ϵ^2`$ (for a disk and planetary system of given dimensions). The imaginary part of equation (47) is a special limit of the conservation law derived in Appendix A.

It is important to note that the self-gravitational term in equation (20) has no effect on the results we have derived. The reason for this is that, if a disk behaves nearly rigidly, the self-precession of the disk is negligibly small even if the mass of the disk exceeds that of the planets.

## 8 Mean-motion resonances

So far we have treated only secular interactions between the planets and disks. However, a planet also interacts with a disk through mean-motion resonances (Lindblad resonances and vertical resonances). The launching of waves at these locations exerts a torque on the planet. If the planet’s orbit is inclined with respect to the resonant annulus of the disk, this inclination will evolve in time.

Borderies, Goldreich, & Tremaine (1984, hereafter BGT) calculated the rate of change of inclination of a satellite due to Lindblad and vertical resonances. They showed that the inclination typically grows exponentially at a rate that depends on the strengths of the resonances present. They assumed implicitly that the launched waves are damped before reaching an edge of the disk from which they might otherwise reflect. In the case of a gas disk, such waves can damp by means of viscous dissipation, radiative damping (Cassen & Woolum 1996), or wave channeling (Ogilvie & Lubow 1999). For narrow gaseous rings between planets, wave damping may well not take place. Instead, reflection from the ring edges may occur, which could greatly reduce the effects of the mean-motion resonances. For simplicity, we will also assume that the waves are damped very close to the resonances and transfer their angular momentum to the disk there.

To incorporate the effect of mean-motion resonances, we modify our equations as follows. (For notational simplicity we here consider only a single planet and a single disk.) The formulae of BGT can be interpreted as giving the horizontal component of the resonant torque between the planet and the disk. We use this to add terms to the planet’s angular momentum equation, in the form

$$m_\mathrm{p}r_\mathrm{p}^2\mathrm{\Omega }_\mathrm{p}\frac{dW_\mathrm{p}}{dt}=\mathrm{}+2\pi \underset{j}{}\frac{Gm_\mathrm{p}^2}{M_{}}s_{\mathrm{r}j}\mathrm{\Sigma }r(W_\mathrm{p}W)|_{r=r_{\mathrm{r}j}},$$
(48)
where $`r_{\mathrm{r}j}`$ is the radius of the $`j`$th resonance, and $`s_{\mathrm{r}j}`$ is the dimensionless strength of the resonance. Similarly, we modify the disk’s angular momentum equation to

$$\mathrm{\Sigma }r^2\mathrm{\Omega }\frac{W}{t}=\mathrm{}+\underset{j}{}\frac{Gm_\mathrm{p}^2}{M_{}}s_{\mathrm{r}j}\mathrm{\Sigma }(WW_\mathrm{p})\delta (rr_{\mathrm{r}j}).$$
(49)
As noted above, this assumes that the waves transfer their angular momentum to the disk very close to the resonances. Provided that the disk is not very cold, the precise location of the wave damping is unimportant because the torque is then communicated through the disk over a large radial extent.

The radii and strengths of the resonances are as follows (BGT); here we have neglected the shifts of the resonant radii caused by precession. Inner vertical resonance ($`m2`$):

$$x=\frac{r_\mathrm{r}}{r_\mathrm{p}}=\left(\frac{m1}{m+1}\right)^{2/3},s_\mathrm{r}=\frac{\pi x^4}{24(m1)}\left[b_{3/2}^{(m)}(x)\right]^2.$$
(50)
Inner Lindblad resonance ($`m2`$):

$$x=\frac{r_\mathrm{r}}{r_\mathrm{p}}=\left(\frac{m1}{m}\right)^{2/3},s_\mathrm{r}=\frac{m\pi x^2}{12(m1)}\left[\left(2m+x\frac{d}{dx}\right)b_{1/2}^{(m)}(x)\right]^2.$$
(51)
Outer vertical resonance ($`m2`$):

$$x=\frac{r_\mathrm{p}}{r_\mathrm{r}}=\left(\frac{m1}{m+1}\right)^{2/3},s_\mathrm{r}=\frac{\pi x^2}{24(m+1)}\left[b_{3/2}^{(m)}(x)\right]^2.$$
(52)
Outer Lindblad resonance ($`m1`$):

$$x=\frac{r_\mathrm{p}}{r_\mathrm{r}}=\left(\frac{m}{m+1}\right)^{2/3},s_\mathrm{r}=\frac{m\pi }{12(m+1)}\left[\left(2m+1+x\frac{d}{dx}\right)b_{1/2}^{(m)}(x)\right]^2.$$
(53)
These quantities are listed in Tables 1 and 2 for all resonances satisfying $`|r_\mathrm{r}r_\mathrm{p}|>0.1r_\mathrm{p}`$. Note that the inner Lindblad resonances act in the opposite sense to the other resonances.

Although equation (49) formally requires the disk tilt $`W`$ to have a cusp at each resonance, these cusps will be very weak in a disk with good radial communication. For a disk that responds nearly rigidly, the effect of the mean-motion resonances is to add a small imaginary part to the planet-disk coupling coefficient $`C^{\mathrm{pd}}`$:

$$C^{\mathrm{pd}}C^{\mathrm{pd}}+2\pi i\underset{j}{}\frac{Gm_\mathrm{p}^2}{M_{}}s_{\mathrm{r}j}\mathrm{\Sigma }_{\mathrm{r}j}r_{\mathrm{r}j}.$$
(54)
The sum over resonances is always positive for an exterior disk but usually negative for an interior disk. According to equation (31), positive sums of this kind tend to cause growth of the non-trivial normal modes. This must compete, however, with the viscous decay implied by equation (47).

## 9 Model for disk and numerical method

We adopt a simplified model for the solar nebula. The disk model is similar to that described in Lubow & Ogilvie (2000). It is characterized by an angular semi-thickness $`H/r=ϵ=\mathrm{constant}`$, and a surface density $`\mathrm{\Sigma }=\mathrm{\Sigma }_0fr^{3/2}`$, where $`\mathrm{\Sigma }_0`$ is a constant and $`f(r)`$ a function that is close to unity in most of the disk, but tapers linearly to zero at the edges (including gaps around planets). The width of the tapers is equal to the local value of $`H`$. The vertical structure is that of a polytrope of index $`n`$, which gives $`=\mathrm{\Sigma }H^2/(2n+3)`$ (as defined in eq. ) if $`H`$ is the true semi-thickness of the polytrope. The basic parameters of the model are then $`ϵ`$, $`\alpha `$, and $`n`$, together with the dimensions and mass of the disk. We normalize the surface density by quoting a nominal disk mass,

$$m_{\mathrm{disk}}=_0^{r_{\mathrm{out}}}\mathrm{\Sigma }_0r^{3/2}\mathrm{\hspace{0.17em}2}\pi r𝑑r=4\pi \mathrm{\Sigma }_0r_{\mathrm{out}}^{1/2},$$
(55)
where $`r_{\mathrm{out}}`$ is the outermost radius. The true disk mass is slightly less than $`m_{\mathrm{disk}}`$ owing to the various edges.

We consider the nebula in the presence of a Jupiter-like planet. As standard parameters, we adopt $`m_\mathrm{p}=0.001M_{}`$, $`m_{\mathrm{disk}}=0.01M_{}`$, $`ϵ=0.075`$, $`\alpha =0.005`$, and $`n=3/2`$. The planet creates a gap in the nebula, leading to two separate disks whose parameters are taken to be $`a_1=0.1r_\mathrm{p}`$, $`b_1=(1g)r_\mathrm{p}`$, $`a_2=(1+g)r_\mathrm{p}`$, $`b_2=20r_\mathrm{p}`$. The standard relative half-width of the gap is taken to be $`g=0.2`$. This situation may be compared with Jupiter in a solar nebula extending out to about 100 AU. The effect of varying the important parameters will be considered below.

Numerical solutions to the equations of Section 5 are obtained by discretizing the disk, converting the integrals into sums, and converting derivatives into centered finite differences, with $`W_{\mathrm{d}k}`$ defined on a set of $`n_\mathrm{r}`$ rings for each disk, and $`𝒢_k`$ on the boundaries between neighboring rings. Planet-disk and disk-disk interactions are computed with softened kernels $`K_m`$, as defined by equation (12). The smoothing length $`h=\mathrm{max}(H_1,H_2)`$ is taken to be the maximum of the semi-thicknesses of the two rings concerned. Equation (49) is implemented by identifying the ring in which each resonance falls, and representing the delta function as $`1/\delta r`$, where $`\delta r`$ is the width of the ring. The problem then reduces to a generalized eigenvalue problem involving a non-Hermitian matrix of dimension $`n_\mathrm{p}+n_\mathrm{d}(2n_\mathrm{r}1)`$. This is solved numerically using the Fortran Sun Performance Library routine ZGEGV on a Sun Ultra 10, using $`n_\mathrm{r}=200`$. Numerical results are quoted in units such that $`G=M_{}=r_\mathrm{p}=1`$. The modes are normalized such that $`W_\mathrm{p}=1`$.

## 10 Numerical results for Jupiter

We first consider a disk with standard parameters, except that it extends only to $`b_2=2`$ instead of $`b_2=20`$. \[Accordingly, we set $`m_{\mathrm{disk}}=(0.1)^{1/2}0.01`$ in this case to fix the same normalization of surface density as in the standard model. The true mass of the disk is only $`0.00179`$.\] This case is easier to understand because the disk is able to maintain good radial communication.

Only two of the non-trivial modes are of interest. The first mode (‘mode 1’) has a precession frequency of $`7.02\times 10^4`$. The inner disk has a nearly rigid tilt $`W20`$ and the outer disk $`W+20`$. If the mean-motion resonances are neglected, this mode has a viscous decay rate of $`4.0\times 10^7`$; when they are included, the mode acquires a net growth rate of $`1.0\times 10^6`$. \[This mode corresponds to equation (36) in rigid disk tilt theory.\]

The second mode (‘mode 2’) has a precession frequency of $`1.46\times 10^3`$. The inner disk has a nearly rigid tilt $`W0.64`$ and the outer disk $`W0.75`$. If the mean-motion resonances are neglected, this mode has a viscous decay rate of $`7.1\times 10^7`$; when they are included, the mode acquires a net growth rate of $`6.6\times 10^6`$.

The remaining non-trivial modes all have significantly larger precession rates and damping rates (with or without mean-motion resonances). They are the proper bending modes of the disks, modified by gravitational coupling to the other components of the system.

When we proceed to the standard model, the numerical solution becomes significantly more complicated. The outer disk extending to $`b_2=20`$ is not able to maintain good radial communication for typical precession rates. The precession frequencies estimated from the nearly rigid theory are comparable to the frequencies of global bending modes in the outer disk. As a result, the nearly rigid theory is no longer accurate for the outer disk.

To simplify the presentation of results, we focus on three modes that are of greatest interest in that they involve the least warping, and have the smallest precession frequencies and damping rates. ‘Mode 1’ and ‘mode 2’ are related to the nearly rigid modes mentioned above, while ‘mode 3’ derives from a bending mode of the outer disk that interferes with the other two when the outer disk fails to behave rigidly. The structures of the three modes are plotted in Fig. 1. The modes include the effects of mean-motion resonances.

In Figs 2 and 3 we plot the precession rates and decay rates of these three modes for the standard model, as we vary each of the parameters $`\alpha `$, $`ϵ`$, $`m_{\mathrm{disk}}`$, and $`g`$ about its standard value. Notice that the variation of decay rate with some parameters is non-monotonic. In particular, the decay rate of mode 1 peaks at $`m_{\mathrm{disk}}0.02`$.

In Appendix A, we derive an expression for the local decay rate of a warped disk due to dissipation. In Fig. 4 we plot for mode 1 the decay rate, as given by equation (A3). Notice that the outer disk dissipation dominates over the inner disk dissipation. Although the dissipation peaks near the planet, it is broadly distributed throughout the outer disk. In the outer disk, the warp $`|rdW/dr|`$ is greatest near the disk radial midpoint and is broadly distributed. Quantity $`|rdW/dr|`$ must vanish at the inner and outer disk edges, due to boundary condition (22). The reason that the peak in Fig. 4 in the outer disk occurs near the planet is that the magnitudes of various base-state quantities, such as $`\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$, are largest there.

In summary, an inclined planet gives rise to large-scale warps. The damping effects of secular interactions dominate over tilt excitation by mean-motion resonances in our model system, provided $`\alpha 0.001`$. For the standard model, the decay time of the longest lived mode, mode 1, is about $`2.7\times 10^5`$ yr. Although this time-scale is short compared to typical disk lifetimes of several million years, it is comparable to or longer than the expected planetary migration time-scale at this radius (Lin, Bodenheimer, & Richardson 1996; Ward 1997).

## 11 Simple estimates

We compare the numerical results with the nearly rigid theory of Section 7 for the standard disk model. To do this, we carry out some rough estimates and drop all factors of order unity. Consistent with our standard model, we assume below that the mass of the planet is less than, or comparable to, the masses of the two disks, and also that the angular momentum of the outer disk is greater than, or comparable to, the angular momenta of the planet and inner disk.

The coupling coefficients are dominated by the parts of the disk closest to the planet, where $`|rr_\mathrm{p}|/r_\mathrm{p}g`$. Here the kernel may be estimated as $`K1/(g^2r)`$. Thus

$$C^{\mathrm{pd}}\frac{1}{g}Gm_\mathrm{p}\mathrm{\Sigma }r$$
(56)
and

$$C^{\mathrm{dd}}G\mathrm{\Sigma }^2r^3.$$
(57)
Here $`\mathrm{\Sigma }`$ is a typical surface density near the planet. The precession rate of a mode can be estimated from either equation (36) or equation (38) as

$$\omega \frac{G\mathrm{\Sigma }}{gr\mathrm{\Omega }}.$$
(58)
We then estimate $`𝒢`$ from equation (20) as

$$𝒢i\omega \mathrm{\Sigma }r^4\mathrm{\Omega }W.$$
(59)
We apply these estimates to equation (47) to obtain an estimate of the viscous damping rate,

$$\mathrm{Im}(\omega )\frac{\alpha }{ϵ^2}\frac{1}{g^2}\frac{(\mathrm{\Sigma }r^2)^2}{M_{}^2}\mathrm{\Omega }.$$
(60)
When mean-motion resonances are taken into account, the planet-disk coupling coefficient acquires a small imaginary part, as in equation (54). The sum over resonances can be estimated using arguments adapted from BGT. For $`g1`$, the number of resonances in the disk scales as $`m1/g`$. The strengths of the Lindblad resonances, which dominate here, scale as $`s_\mathrm{r}m^2`$. Thus

$$\mathrm{Im}(C^{\mathrm{pd}})\frac{1}{g^3}\frac{Gm_\mathrm{p}^2\mathrm{\Sigma }r}{M_{}}.$$
(61)
¿From equation (36) this would provide a resonant growth rate

$$\mathrm{Im}(\omega )\frac{1}{g^3}\frac{m_\mathrm{p}\mathrm{\Sigma }r^2}{M_{}^2}\mathrm{\Omega }.$$
(62)
The ratio of viscous decay to resonant growth can then be estimated as

$$\frac{\mathrm{decay}}{\mathrm{growth}}\left(\frac{\alpha }{ϵ^2}\right)\left(\frac{g\mathrm{\Sigma }r^2}{m_\mathrm{p}}\right).$$
(63)
The first factor depends on the viscosity and temperature of the disk. The second is related to the ratio of the disk mass within (say) two gap widths of the planet to the mass of the planet itself. For our standard parameters, these first factor is close to unity, while the second is somewhat less than unity. This favors tilt growth, provided that the disk indeed behaves nearly rigidly.

## 12 Connection with WKB bending wave theory

Disks subject to forcing by a misaligned companion have been studied in the context of planetary rings, but with strikingly different results. Within Saturn’s A ring, tightly wrapped waves are launched from a mean-motion vertical resonance by a misaligned satellite, as described by Shu, Cuzzi, & Lissauer (1983). Furthermore, the torque carried by the waves is independent of $`\alpha `$ and $`ϵ`$, unlike the present case.

As emphasized earlier, the difference lies in the fact that the waves here are of low frequency such that the parameter $`\omega r/c_\mathrm{s}`$ is typically less than unity for protostellar disks with secular resonances. For the mean-motion vertical resonances in planetary rings, this quantity is much greater than unity.

Nonetheless, a formal connection between the equations used in this paper and those used for planetary rings can be made by considering an artificially cold protostellar disk, such that $`\omega r/c_\mathrm{s}`$ is much greater than unity. Furthermore, we consider the case that the disk mass is much greater than the planetary mass. Appendix B outlines the derivation of the dynamical equations in the WKB limit. In that limit we obtain

$$D(r)W2\pi i\sigma _kG\mathrm{\Sigma }\frac{dW}{dr}+\left(\frac{}{2\mathrm{\Sigma }}\right)\frac{\mathrm{\Omega }^3}{[\omega +(\kappa ^2\mathrm{\Omega }^2)/(2\mathrm{\Omega })i\alpha \mathrm{\Omega }]}\frac{d^2W}{dr^2}\frac{2Gm_\mathrm{p}}{r^2}K(r_\mathrm{p},r)W_\mathrm{p},$$
(64)
where $`\sigma _k=\pm 1`$ with the sign depending on whether the waves are trailing or leading, and

$$D(r)=2\omega \mathrm{\Omega }+\frac{2Gm_\mathrm{p}}{r^2}K(r_\mathrm{p},r).$$
(65)
For $`\omega \mathrm{\Omega }`$, we can re-express $`D(r)`$ in the standard form

$$D(r)=\mu ^2(\mathrm{\Omega }^{}\omega )^2,$$
(66)
where

$$\mu ^2=\frac{^2\mathrm{\Phi }^{}}{z^2}|_{z=0},r\mathrm{\Omega }^2=\frac{\mathrm{\Phi }^{}}{r}|_{z=0}.$$
(67)
Potential $`\mathrm{\Phi }^{}`$ is due to the star and axisymmetric contributions of the planet, but does not include contributions due to the disk. This equation is identical in form to equation (20) of Shu et al. (1983) for $`m=1`$, apart from the pressure term, the third term in the equation (64).

The pressure term arises because the warping of the disk introduces horizontal pressure gradients. The resulting shearing horizontal motions give rise to the hydrodynamic torque $`2\pi 𝑮`$, which is the dominant factor in the propagation of bending waves in protostellar disks. The warping is described by the $`n=0`$ mode in the notation Lubow and Pringle (1993). The dispersion relation for adiabatic oscillations of this mode in an exactly Keplerian disk, whose unperturbed state is vertically isothermal, follows from equation (54) of Lubow and Pringle (1993) with dimensionless frequencies $`F=(\mathrm{\Omega }\omega )/\mathrm{\Omega }1`$ and $`\kappa =1`$. In the low-frequency regime ($`|\omega |\mathrm{\Omega }`$), one obtains then that $`\omega =\pm kH\mathrm{\Omega }/2`$, independent of the exponent for adiabatic compression $`\gamma `$. Although the form of the dispersion relation is reminiscent of that for a compressive acoustic mode, this mode is largely incompressible, as suggested by the lack of dependence on $`\gamma `$ in the dispersion relation. The same dispersion relation (as was obtained by Papaloizou & Lin 1995) follows from equation (64), with $`W`$ varying as $`\mathrm{exp}(ikr)`$, and $`=(c_\mathrm{s}/\mathrm{\Omega })^2\mathrm{\Sigma }`$, and by taking $`G=0`$ (ignoring gravity), $`\kappa =\mathrm{\Omega }`$, and $`\alpha =0`$.

For this limit where the planet mass is small compared to the disk mass, the resonance condition $`D(r)=0`$ becomes

$$\frac{m_\mathrm{p}}{m_\mathrm{d}}\frac{K(r,r_\mathrm{p})}{K}\left(\frac{r_\mathrm{p}}{r}\right)^{1/2}=1$$
(68)
where $`K`$ denotes the mass weighted average of $`K`$ over the disk (see Appendix B for more details). For the standard disk model, the resonance condition is satisfied at radial distances of about 0.7 and 1.3 within the inner and outer disks respectively.

This system contains a secular particle resonance where the precession rate of a particle matches the precession rate of the disk-planet system. A free particle (e.g. an asteroid) will precess due to the gravitational effects of both the disk and the planet. That is, a free particle resonance will occur when $`D(r)=0`$ in equation (66), but $`\mathrm{\Phi }^{}`$ now has an additional contribution due to the disk. The particle resonance location is therefore not the same as the disk resonance location. In other words, in locating a particle secular resonance, the gravity of the gas disk must be included. Particles oscillating vertically in a non-oscillating gas disk have an additional vertical restoring force that alters their precession frequency and shifts their resonance site. This is a special property of the $`m=1`$ bending modes (tilt modes), which is a consequence of the fact that there is no self-precession in a rigidly tilted disk. Near resonance, a tilt mode essentially behaves rigidly, due to its locally long wavelength.

In the case of a planetesimal (particle) disk, the nature of the pressure and viscous interaction may well differ from that of a gaseous disk. On the other hand, the secular resonance condition given by equation (68) likely applies to a low-mass planet interacting with a cold planetesimal disk.

## 13 Planet in central hole

Current models for the early evolution of close-orbiting extra-solar planets, such as 51 Peg b (Mayor & Queloz 1995) with an orbital radius of about 0.05 AU, involve the planet migrating into a central hole in the disk (Lin et al. 1996). The planet is envisioned to migrate into the hole by disk tidal forces until its 2:1 outer Lindblad resonance lies just inside the inner edge of the disk. The planet would then remain at this radius, without further migration.

We consider the tilt stability involving both secular interactions and mean-motion resonances of such a planet-disk system, based on the standard model parameters in Section 9. In this case, there is only an outer disk of mass 0.01 $`M_{}`$, which extends from just outside the planet’s 2:1 outer Lindblad resonance $`1.31r_\mathrm{p}`$ (0.066 AU) to outer radius 2000 $`r_\mathrm{p}`$ (100 AU). Other disk parameters are the same as in the standard model ($`ϵ=0.075`$, $`\alpha =0.005`$, and $`n=3/2`$). The only mean-motion resonance present in the disk is the rather weak 1:3 resonance (2.08 $`r_\mathrm{p}`$), which acts to increase inclination.

Fig. 5 shows the eigenfunction of the lowest mode of the system. The eigenfrequency is $`5.11\times 10^7+4.6\times 10^8i`$ in units of $`\mathrm{\Omega }_\mathrm{p}`$, corresponding to a decay time-scale of only $`3.9\times 10^4`$ yr. The inner part of the disk interacts nearly rigidly with the planet, owing to good radial communication $`\omega r_\mathrm{p}/c_\mathrm{s}1`$ (see lower panel of Fig. 5). However, the outer part of the disk is not in good communication ($`\omega r/c_\mathrm{s}`$ is of order unity), resulting in a warp and dissipation.

## 14 Disk and two planets

We consider here the case of Jupiter and Saturn interacting with the portion of the solar nebula that lies interior to the orbit of Jupiter, using the equations of Section 5. This case is of interest because the classical $`\nu _{16}`$ nodal secular resonance lies within the disk. We ignore the effects of mean-motion resonances.

The rigid tilt model involving three objects (here the two planets plus disk) was considered at the end of Section 6. As is mentioned there, a secular resonance of the two planets with the disk is possible, if the disk is low in mass compared with the planets. If the disk mass is greater than the planet masses, then the planets individually interact with the disk in a manner similar to the case in Section 10.

We have computed the modes of the system numerically. We consider Jupiter (planet 1) and Saturn (planet 2) in circular orbits with the present values of mass ($`m_{\mathrm{p1}}=0.0009545`$, $`m_{\mathrm{p2}}=0.0002858`$) and semi-major axis ($`r_{\mathrm{p2}}=1.833`$). A test particle would experience nodal secular resonances at semi-major axes $`0.381`$ and $`2.39`$, corresponding to $`1.98`$ AU and $`12.4`$ AU (Murray & Dermott 1999). We consider a partially depleted solar nebula interior to Jupiter’s orbit. Our standard model has $`m_{\mathrm{disk}}=0.001`$, $`ϵ=0.075`$, $`\alpha =0.005`$, $`n=3/2`$, and disk inner and outer radii $`a_1=0.1`$ and $`b_1=1g=0.8`$, respectively, in units of Jupiter’s orbital radius. The inner secular resonance for a test particle therefore lies inside the disk.

For these parameters, the numerically determined eigenfrequencies of the two lowest modes are $`2.27\times 10^4+6.0\times 10^9i`$ and $`7.66\times 10^4+5.1\times 10^7i`$. The response of the disk is very nearly rigid and nothing special happens at the location of the secular particle resonance. This is to be expected, since $`\omega r/c_\mathrm{s}1`$ throughout the disk. The effect of self-gravitation is negligible even though the mass of the disk is comparable to that of the planets.

If a low-mass disk is extremely cold, it is unable to establish radial communication very effectively on secular time-scales. It cannot respond nearly rigidly, but instead launches a wave at the location of the particle resonance. (The rigid response is essentially the limit of a wavelike response when the wavelength becomes large compared to the size of the disk.) In Fig. 6 we show the wavelike response generated at the $`\nu _{16}`$ secular resonance at $`2`$ AU. A very thin disk ($`ϵ=10^5`$) is required to see this effect, and even then the wavelength is not much shorter than the size of the disk (so a WKB or tight-winding treatment of the wave might not be particularly accurate). We set $`\alpha =ϵ=10^5`$ so that the wave would be partially damped on reaching the inner boundary. We also chose a low disk mass, $`m_{\mathrm{disk}}=10^6`$, so that the disk has a Toomre parameter $`Q>1`$. The eigenfrequency for this lowest mode is $`2.35\times 10^4+1.6\times 10^8i`$, corresponding to a damping time of about $`1.1\times 10^8`$ yr. The resonantly launched waves damp inclination faster than those in the warmer disk model, even though the disk mass is much lower.

## 15 Summary

We have carried out a normal mode analysis of planet-disk systems that includes the effects of secular interactions and mean-motion resonances. The planets and disks interact gravitationally, and the disks communicate internally through gas pressure, self-gravity, and turbulent viscosity. Secular interactions of misaligned planet-disk systems give rise to global effects in the disks, owing mainly to the effects of gas pressure. The low frequencies $`\omega `$ associated with secular modes allow pressure to communicate effectively over radial distances in disks that are comparable to (and somewhat greater than) the planet’s orbital radius. Over such distances, disks behave rigidly with little dissipation. On the other hand, over larger distances ($`c_\mathrm{s}/\omega `$), which may arise in a continuous disk exterior to the planet’s orbit, the radial communication breaks down and large-scale warps occur, along with enhanced dissipation (see Figs 1 and 5). The protostellar disk response to secular interactions typically falls between that of a global nearly rigidly tilted disk (see Section 7) and that of a locally launched wave (see Sections 12 and 14).

Secular interactions act to decrease inclination by means of turbulent dissipation of horizontal shearing motions within the warp for a simple $`\alpha `$ disk prescription. The dissipation of the warp is peaked somewhat near the planet, but is broadly distributed throughout the disk (Fig 4). Secular interactions are in competition with the effects of mean-motion resonances (BGT), which can act to increase inclination. For standard disk parameters, the secular interactions involving a single planet generally dominate and the inclination decays in time. For a young Jupiter interacting with the solar nebula, tilt decay takes places for $`\alpha 0.001`$, with a typical alignment time-scale of order $`10^5`$ yr (see Fig. 3). For a close-orbiting planet in the central hole of the disk, the alignment is more rapid. These results suggest that a planet formed within a disk will remain coplanar with the disk, at least as a consequence of disk interactions.

The classical $`\nu _{16}`$ particle secular resonance, due to Jupiter and Saturn, is expected to have resided within the solar nebula. Again, due to pressure effects, this resonance is greatly broadened in the gas disk. Consequently, the nebula responds nearly rigidly, with some decay in inclination. A very cold planetesimal disk might exhibit a genuinely wavelike response (cf. Fig 6).

The analysis in this paper has focused mainly on the effects of a single planet interacting with its interior and exterior disks. For a system of a multiple planets that opens multiple disk gaps, the situation is more complicated because of reflections that occur at the disk edges associated with the gaps. Such gaps can decrease the effective size of the disk region over which a planet can interact, thereby resulting in a weaker warp and weaker dissipation. The hydrodynamic torque cannot propagate across a gap, while the gravitational torque may at least partially do so. It is quite possible that much longer lived noncoplanar modes may result if the outer part of the nebula is almost disconnected in this way.

This work was supported by NASA grants NAG5-4310 and NAG5-10732, the STScI visitor program, and the Institute of Astronomy visitor program. GIO acknowledges the support of Clare College, Cambridge through a research fellowship and the Royal Society through a University Research Fellowship.

## Appendix A Derivation of the local decay rate of bending disturbances

Starting from equations (13)–(15) we obtain the conservation law

$$\frac{d}{dt}(L_z)=\underset{k}{}_{a_k}^{b_k}\frac{4\alpha |𝒢_k|^2}{_kr^4\mathrm{\Omega }^2}2\pi r𝑑r,$$
(A1)
where

$$L_z=\underset{i}{}\frac{1}{2}m_{\mathrm{p}i}r_{\mathrm{p}i}^2\mathrm{\Omega }_{\mathrm{p}i}|W_{\mathrm{p}i}|^2+\underset{k}{}_{a_k}^{b_k}\left(\frac{1}{2}\mathrm{\Sigma }_kr^2\mathrm{\Omega }|W_{\mathrm{d}k}|^2+\frac{2|𝒢_k|^2}{_kr^4\mathrm{\Omega }^3}\right)2\pi r𝑑r.$$
(A2)
The quantity $`L_z`$ is the vertical angular momentum associated with the bending disturbance. When an initially horizontal ring of matter having orbital angular momentum $`J`$ is tilted through an angle $`\beta =|W|`$, the change in the vertical component of angular momentum is $`J(\mathrm{cos}\beta 1)\frac{1}{2}J|W|^2`$ when $`|W|1`$. The term proportional to $`|𝒢|^2`$ represents an additional angular momentum stored in the bending wave. For a short-wavelength (WKB) bending wave in a non-self-gravitating disk, both terms in the integral for $`L_z`$ are equal. For a nearly rigidly tilted disc, the first term in the integral dominates. In the appropriate limit, this first term agrees with the expression for the angular momentum of a bending wave given by Bertin & Mark (1980, eq. \[C12\]).

In an inviscid system $`L_z`$ is conserved. When viscosity is present, the positive definite quantity $`L_z`$ decays monotonically to zero. The negative angular momentum of the bending disturbance is transferred to the disk through viscosity, and causes accretion. (It is nevertheless consistent that we neglect the time-dependence of $`\mathrm{\Sigma }`$ in a linear theory, because the induced accretion rate is quadratic in $`W`$.)

We now define a local damping rate $`\gamma (r)`$ by

$$\gamma =\left(\frac{1}{L_z}\right)\frac{2\alpha |𝒢|^2}{r^4\mathrm{\Omega }^2}2\pi r\underset{k}{}(b_ka_k).$$
(A3)
This has the property that

$$\frac{\underset{k}{}_{a_k}^{b_k}\gamma 𝑑r}{_k(b_ka_k)}=\frac{1}{2}\frac{d}{dt}\mathrm{ln}(L_z).$$
(A4)
In the case of a normal mode, the radial average of $`\gamma `$ corresponds to the damping rate $`\mathrm{Im}(\omega )`$ of the mode.

## Appendix B Derivation of the bending wave equations in the WKB limit

In this Appendix we derive equations (64) and (68) of the text. We consider equations (19)–(21) for a single planet interacting with a single disk.

$`i\omega m_\mathrm{p}r_\mathrm{p}^2\mathrm{\Omega }_\mathrm{p}W_\mathrm{p}`$ $`=`$ $`{\displaystyle Gm_\mathrm{p}\mathrm{\Sigma }K(r,r_\mathrm{p})i(WW_\mathrm{p})\mathrm{\hspace{0.17em}2}\pi r𝑑r},`$ (B1)
$`i\omega \mathrm{\Sigma }r^2\mathrm{\Omega }W`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{d𝒢}{dr}}+Gm_\mathrm{p}\mathrm{\Sigma }K(r,r_\mathrm{p})i(W_\mathrm{p}W)`$ (B2)
$`+{\displaystyle G\mathrm{\Sigma }\mathrm{\Sigma }^{}K(r,r^{})i(W^{}W)\mathrm{\hspace{0.17em}2}\pi r^{}𝑑r^{}},`$
$$i\omega 𝒢+\left(\frac{\kappa ^2\mathrm{\Omega }^2}{\mathrm{\Omega }^2}\right)\frac{i\mathrm{\Omega }}{2}𝒢+\alpha \mathrm{\Omega }𝒢=\frac{r^3\mathrm{\Omega }^3}{4}\frac{dW}{dr}.$$
(B3)
We now apply the WKB approximation for a cold disk, $`\omega rc_\mathrm{s}`$. In this regime we assume

$$H^1\left|\frac{\mathrm{ln}W}{r}\right|r^1.$$
(B4)
The requirement on $`H^1`$ is used in the derivation of the warp equations (see Papaloizou & Lin 1995; Lubow and Ogilvie 2000). From equation (B3), we obtain

$$\frac{1}{r}\frac{d𝒢}{dr}\frac{i}{4}\frac{r^2\mathrm{\Omega }^3}{[\omega +(\kappa ^2\mathrm{\Omega }^2)/(2\mathrm{\Omega })i\alpha \mathrm{\Omega }]}\frac{d^2W}{dr^2}.$$
(B5)
For the self-gravity term, we have that

$$G\mathrm{\Sigma }\mathrm{\Sigma }^{}K(r,r^{})i(W^{}W)\mathrm{\hspace{0.17em}2}\pi r^{}𝑑r^{}\pm \pi G\mathrm{\Sigma }^2r^2\frac{dW}{dr}.$$
(B6)
To obtain equation (B6), we recognize that self-gravity is dominated by local contributions in the WKB limit and apply the approximation that $`K(r,r^{})r/[2\pi (r^{}r)^2]`$ for $`r^{}r`$. The integral is computed by extending $`r^{}`$ to the complex plane and integrating over a closed contour. The contour begins along the negative real axis, includes a small semicircle around the singularity at $`r^{}=r`$, continues along the real axis, and closes on itself through a large semicircle. The large semicircle resides in either the upper or lower half-plane, depending on whether the waves are trailing or leading. Equation (B6) then follows from the residue theorem. This choice of trailing or leading waves determines the sign of the result. Substituting equations (B5) and (B6) into equation (B2), we obtain equation (64).

Integrating equation (B2) over the disk area and using equation (B1), we have that

$$J_\mathrm{d}W=J_\mathrm{p}W_\mathrm{p},$$
(B7)
where $`W`$ denotes the angular momentum weighted average of $`W`$. In the limit that $`J_\mathrm{p}J_\mathrm{d}`$, we expect that $`W`$ on the right hand side of equation (B1) can be ignored after integration. It then follows from equation (B1) that

$$\omega =\frac{G\mathrm{\Sigma }K(r,r_\mathrm{p})2\pi r𝑑r}{r_\mathrm{p}^2\mathrm{\Omega }_\mathrm{p}}=\frac{C^{\mathrm{pd}}}{J_\mathrm{p}}.$$
(B8)
Applying this equation for $`\omega `$ to the resonance condition $`D(r)=0`$, for $`D`$ defined by equation (65), we obtain the resonance condition (68).


text:
# An Empirical Ultraviolet Iron Spectrum Template Applicable to Active Galaxies

## 1. Motivation

This work was motivated by our desire to study the intrinsically emitted ultraviolet (UV) broad line profiles of active galactic nuclei (AGNs) to get more accurate and reliable information on e.g., the chemical abundances and dynamical and ionization structure in their broad line region.

In addition to the commonly studied strong, permitted lines of Ly $`\alpha \lambda `$ 1216, Si iv+O iv\] $`\lambda `$ 1400, C iv$`\lambda `$ 1549, C iii\] $`\lambda `$ 1909, and Mg ii$`\lambda `$ 2800, AGN restframe-UV spectra also contain weaker lines such as various transitions of carbon, silicon, oxygen, nitrogen, helium, and iron. The importance of these lines lies not only in their presence and strengths but also in their effective “contamination” of the stronger lines. Iron is by far the strongest of these “contaminants” of the UV (and optical) emission lines. Due to the large number of electron levels in iron atoms, thousands of emission-line transitions of several ionization states are distributed throughout the UV and optical spectral regions (e.g., Kurucz 1997). The weak lines of iron and other metals blend together, especially when the AGN intrinsic line width is large. This heavy blending of weak, broad lines, dominated by iron transitions, forms a pseudo-continuum above the intrinsically emitted continuum level (Wills, Netzer & Wills 1985) even in high signal-to-noise AGN spectra. The presence of such a pseudo-continuum severely complicates the study of the wings of the prominent resonance lines, the weak line features themselves, and the intrinsically emitted continuum emission. It is important to make reliable measurements of these features in order to study issues such as continuum – line relationships, dynamical and ionization structure, and chemical abundances in the AGN broad line region.

Our preliminary analysis indicates that line measurement uncertainties of order $``$5 – 25% are expected in the prominent lines if the iron emission is not corrected for. Some weak emission lines may, however, be so severely blended with iron emission that the uncertainty approaches 100%.

## 2. Why an Empirical Template?

A viable way to correct for the contaminating iron emission is to generate an empirical iron emission template which can be used to ‘model’ the iron emission in individual AGN spectra. An alternative is to theoretically model the iron emission. However, to date such modeling has only been moderately successful mainly due to inadequate atomic data and limited understanding of the rather complicated iron emission mechanism. Much improved atomic data have become available in recent years (e.g., Hummer et al. 1993; Seaton et al. 1994; Kurucz 1997) allowing much more advanced theoretical models to be developed. This work is currently in progress (e.g. E. Verner 2001, private communication; Verner et al. 1999). An empirical iron template is thus extremely useful until adequate theoretical models and an understanding of the iron spectrum emitted by AGNs are available. A template is not model dependent but comprises the iron emission transitions and line ratios as emitted by an active galaxy. It permits the fitting and subtraction of the iron emission in individual AGN spectra, as the template can be broadened and scaled to match the iron emission in each AGN spectrum. In addition, this allows us to create a database of (fitted) iron spectra which can help constrain theoretical models of AGN iron emission, thereby helping us to better understand this emission. An optical template, based on data of I Zw 1, has already been successfully applied to the bright quasar survey (Boroson & Green 1992).

## 3. The Template

We present a UV Fe ii and Fe iii empirical template ranging from 1250 – 3090 Å based on high quality HST archival spectra of the nearby AGN, I Zw 1. The data were originally obtained, presented, and analyzed by Laor et al. (1997), who also identify most of the emission and absorption features in the spectra. The UV spectrum of I Zw 1 was chosen for several reasons. The spectrum is rich and strong in the iron emission, the prominent resonance lines are relatively narrow (full width at half maximum $``$ 900 km s-1) and an optical iron emission template based on this object already exists, as mentioned.

A strong and rich iron spectrum permits the detection of weak iron features and the identification of as many iron features as are presumably typically emitted by AGNs. The relatively narrow lines allow us to deblend the line emission in the process of generating the iron spectrum. In order to isolate the iron emission, the other emission and absorption lines were modeled (composite Gaussian profiles were used to reproduce the line profile shapes) and subtracted.

The narrow I Zw 1 line width reveals that a number of line-emission transitions may dominate the C iii\] $`\lambda `$ 1909 line blend; the C iii\] contribution is small for I Zw 1. Hence, caution must be exercised when studying C iii\] as direct measurements of this line can be significantly in error for broader lined active galaxies where the emission lines blend together and are difficult to discern. This affects general line studies, broad-line region density determinations, and abundance studies. This heavy blending of C iii\] emphasizes the motivation for this work.

The Fe ii and Fe iii emission templates are shown in Figure 1. Small residuals of the subtracted line fits are set to zero to avoid introducing unnecessary uncertainties into the AGN spectra to which the templates are applied. In the region around Ly $`\alpha `$ the iron emission is either too weak to be reliably isolated or non-existent, except in the red wing of Ly $`\alpha `$. Blueward of Ly $`\alpha `$ the Lyman forest heavily complicates a reliable isolation of the iron emission. The strength of the Fe iii $`\lambda `$ 1540Å component is very uncertain as it is heavily blended with C iv emission, and blueshifted C iv emission is possibly present. Vestergaard & Wilkes (2001) present details of the template generation, its application, and address the limitations of the template fitting method.

## 4. Future Work

The Fe ii and Fe iii templates are used to fit the UV iron emission in large samples of quasar spectra (Forster et al. 2001; Vestergaard et al. 2001, in preparation) and can be used to fit any AGN UV spectrum with intrinsic line width $``$900 km s-1.

The fitting and subtraction of the contaminating iron emission in UV spectra of AGN will improve studies of their broad emission lines and will help improve photoionization modeling of the broad line region. The database of fitted iron emission spectra from individual active galaxies can then be used to attempt to understand this iron emission in terms of physical models.

### Acknowledgments.

We are grateful to Dima Verner for discussions and direction to updated iron line lists and to Ari Laor for useful discussions. MV is very pleased to thank the Smithsonian Astrophysical Observatory, where most of this work was done, for their hospitality and gratefully acknowledges financial support from the Columbus Fellowship, the Danish Natural Sciences Research Council (SNF-9300575), the Danish Research Academy (DFA-S930201), and a Research Assistantship at Smithsonian Astrophysical Observatory (NAGW-4266, NAGW-3134, NAG5-4089). BJW gratefully acknowledges financial support from NASA contract NAS8-39073 (Chandra X-ray Center).

## References

Boroson, T. A., & Green, R. F. 1992, ApJS, 80, 109

Forster, K., Green, P.J., Aldcroft, T.A., Vestergaard, M., Foltz, C., & Hewitt, P., 2001, ApJS, 134, 35

Hummer, D. G. et al. 1993, A&A, 279, 298

Kurucz, M. 1997, CD-ROM 23 (Cambridge: SAO)

Laor, A., Jannuzi, B.T., Green, R. F., & Boroson, T. A. 1997, ApJ, 489, 656

Seaton, M. J., Yan, Y., Mihalas, D., & Pradhan, A. K. 1994, MNRAS, 266, 805

Verner, E. M. et al. 1999, ApJS, 120, 101

Vestergaard, M., & Wilkes, B.J. 2001, ApJS, 134, 1

Wills, B.J., Netzer, H., & Wills, D. 1985, ApJ, 288, 94