text:
# Theoretical estimates for the largest Lyapunov exponent of many-particle systems

## I Introduction

The largest Lyapunov exponent measures the sensitivity to initial conditions in a dynamical system. In low-dimensional models, the Lyapunov exponent sets a limit to the prediction of the time evolution of a given state of the system. In the high dimensional systems encountered in thermodynamics, one abandons from the start a description in terms of single (microscopic) states, and resorts to a statistical approach. In these cases a positive Lyapunov exponent is usually welcome, as it is a necessary condition for the validity of the Boltzmann-Gibbs scenario. A Lyapunov exponent that becomes null in the thermodynamic limit is a signal of anomalous behaviors. For instance, metastable phases of some long-range interacting systems, where the Lyapunov exponent vanishes in the large $`N`$ limit, exhibit breakdown of ergodicity, anomalous diffusion, and non-Maxwell velocity distributions . An extension of standard statistical mechanics is required for the theoretical explanation of such phenomena .

In many-particle systems the Lyapunov exponent is also an indicator (order parameter) of phase transitions  and may be related to transport coefficients .

In practice, for a given system, the largest Lyapunov exponent must be obtained through numerical simulations, typically using the method developed by Benettin and collaborators  (a proof that this method gives the Lyapunov exponent of Oseledec’s theorem  can be found in Ref. ). For some special cases, e.g., hard-sphere systems, the theory of Lyapunov exponents is remarkably developed . This contrasts with the situation in smooth Hamiltonian systems, where comprehensive analytical estimates are scarce. In 1984, when studying two-dimensional billiards, Benettin made the first step towards the construction of random matrix methods for modeling universal features of Lyapunov spectra . These methods were further developed by other authors and applied to many-particle smooth Hamiltonians .

To our knowledge, the first theory for estimating the largest Lyapunov exponent of a specific Hamiltonian system, and its dependence on the system parameters, was formulated by Pettini and coworkers some years ago . In this approach, the dynamics is geometrized by absorbing the force terms into a suitable metric, thus mapping the Hamiltonian problem onto a geodesic motion on a curved manifold. After making the “quasi-isotropy” approximation, the Hamilton equations for the tangent vectors become decoupled. As a consequence, the initial system of $`2N`$ differential equations has been reduced to only two equations. While the original problem was governed by the Hessian matrix of the potential, of size $`N\times N`$, the new (reduced) one is controlled by the Laplacian of the potential, $`\mathrm{}𝒱(t)`$, a scalar function of time. Thereafter, $`\mathrm{}𝒱(t)`$ is treated as Gaussian white noise and the $`2\times 2`$ system of differential equations is solved using the methods developed by van Kampen and others . See Ref.  for a review.

When applied to a Fermi-Pasta-Ulam chain , the so-called “geometric method” was extremely successful in reproducing the largest Lyapunov exponent over the entire energy range . However, in other cases the agreement is not so good. For instance, in a chain of rotators with first-neighbor (bounded) interactions, the method works well only in the low and high energy regimes, where the dynamics is weakly chaotic (integrable in the limits $`E0,\mathrm{}`$). In the intermediate region of stronger chaos, the theory has to be amended to obtain a good agreement with simulations . This and other examples  raise several questions concerning the domain of validity of the theory. What is the nature of the quasi-isotropic approximation? Or, what are the parameters that control the quality of the estimates of the theory? Is the geometric method perturbative? If so, what are the next leading corrections?

In this paper we present an alternative theoretical approach in which the validity domains of the successive approximations can be precisely delimited. The basic idea is to employ van Kampen’s methods  to solve the original system of $`2N`$ differential equations for the evolution of tangent vectors. By applying this scheme to a three dimensional dilute gas, Barnett et al.  established a link between the Lyapunov exponent and the self-diffusion coefficient (see also ). We show that this approach can be extended to reach other systems, like Fermi-Pasta-Ulam chains and lattices of classical spins, either with short or long-range interactions. By doing so, we shall settle down a connection with the results of the geometric method and suggest some answers to the above-mentioned questions.

The paper has been organized as follows. Section II presents the theory that leads to an estimate of the largest Lyapunov exponent. This is a perturbative theory, which rests on a cummulant expansion. We argue that the general (perturbative) solution can be obtained by diagonalizing a small-dimension matrix. In Sect. III we analyze an approximation that reduces the problem to diagonalizing a $`3\times 3`$ matrix. Section IV discusses some examples which illustrate the working of the theory. The connection between our results and those obtained by the geometric method is discussed in Sect. V. Finally, Sect. VI contains the concluding remarks.

## II Theory

The theory we present in this section can, in principle, be applied to any smooth Hamiltonian system. For simplicity, and for the sake of comparisons with the geometric method, we restrict ourselves to “natural” Hamiltonians,

$$=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2m}+𝒱(q_1,\mathrm{},q_N),$$
(1)
where $`q_i`$ and $`p_i`$, are conjugate position-momentum coordinates. Other Hamiltonians can be considered but they may require modifications of the theory.

The Hamilton equations can be written in the compact form

$$\dot{x}=𝐉\frac{}{x},$$
(2)
where we have introduced the $`2N`$-dimensional column vector $`x`$, $`x=(q_1,\mathrm{},q_N,p_1,\mathrm{},p_N)^T`$, the superscript meaning “transposed”, and the symplectic matrix $`𝐉`$,

$$𝐉=\left(\begin{array}{cc}0& \\ 𝟙& 0\end{array}\right),$$
(3)
with $`𝟙`$ the $`N\times N`$ identity matrix. Differentiating the Hamilton equations, one obtains the evolution equations for tangent vectors $`\xi =(\delta q_1,\mathrm{},\delta q_N,\delta p_1,\mathrm{},\delta p_N)^T`$,

$$\dot{\xi }=𝐀(t)\xi .$$
(4)
For a Hamiltonian of the special form (1), and setting $`m=1`$, the operator $`𝐀`$ has the simple structure

$$𝐀(t)=\left(\begin{array}{cc}0& \\ 𝐕(t)& 0\end{array}\right).$$
(5)
Here $`𝐕`$ is the Hessian matrix of the potential $`𝒱`$, namely

$$V_{ij}=\frac{^2𝒱}{q_iq_j}.$$
(6)
Once initial conditions $`x_0`$ and $`\xi _0`$ have been specified, Eqs. (2) and (4) allow one to find the Lyapunov exponent $`\lambda `$ by calculating the limit

$$\lambda =\underset{t\mathrm{}}{lim}\frac{1}{2t}\mathrm{ln}|\xi (t;x_0,\xi _0)|^2.$$
(7)
We will assume that for any initial condition $`x_0`$, the phase-space trajectory $`x(t;x_0)`$ is ergodic on its energy shell. This implies that $`\lambda `$ is independent of initial conditions $`x_0`$, which can then be chosen randomly according to the microcanonical distribution. There will also be no dependence on initial tangent vectors, because if $`\xi _0`$ is also chosen randomly, it will have a non-zero component along the most expanding direction.

If the corrections to the exponential law in Eq. (7) go to zero fast enough as $`t\mathrm{}`$, one can also write

$$|\xi (t;x_0,\xi _0)|^2e^{2\lambda t},$$
(8)
brackets meaning microcanonical averages over $`x_0`$. We will prefer the estimate of Eq. (8) because the averaging procedure is crucial for finding an analytical expression for the Lyapunov exponent. In case of doubt, the equality of the exponents defined by Eqs. (7) and (8) can be tested numerically, e.g., using the data generated by Benettin’s algorithm.

By letting $`x_0`$ be a random variable, a stochastic process $`𝐕(t;x_0)`$ is defined, and Eq. (4) can be thought as a stochastic differential equation. However, the quantity we are interested in is the square of the norm of $`\xi `$, which can be written as the trace of the “density matrix” $`\xi \xi ^T`$. Thus, we must focus on the equation for the evolution of $`\xi \xi ^T`$:

$$\frac{\mathrm{d}}{\mathrm{d}t}(\xi \xi ^T)=𝐀\xi \xi ^T+\xi \xi ^T𝐀^T\widehat{𝐀}\xi \xi ^T,$$
(9)
the rightmost identity defining the linear superoperator $`\widehat{𝐀}`$. Except for the fact that we must deal now with a superoperator, Eq. (9) is not different from Eq. (4) and can be handled with the same techniques. For the purpose of the perturbative approximations that will follow, the operator $`\widehat{𝐀}`$ is split into two parts

$$\widehat{𝐀}=\widehat{𝐀}_0+\widehat{𝐀}_1(t),$$
(10)
where $`\widehat{𝐀}_0`$ corresponds to the evolution in the absence of interactions. In our case $`\widehat{𝐀}_0`$ and $`\widehat{𝐀}_1`$ are associated with

$$𝐀_0=\left(\begin{array}{cc}0& \\ 0& 0\end{array}\right)\text{and}𝐀_1=\left(\begin{array}{cc}0& 0\\ 𝐕(t)& 0\end{array}\right),$$
(11)
respectively. Whenever $`𝐀_1(t)`$ is small (in a sense that will be discussed below), it is possible to manipulate Eq. (9) to derive an explicit expression for the evolution of the average of $`\xi \xi ^T`$. A clear exposition of this derivation, together with a very detailed discussion of its domain of validity has been given by van Kampen . We just outline the basic steps: (a) Rewrite Eq. (9) in the interaction representation associated with $`\widehat{𝐀}_0`$. (b) Write the propagator as a time ordered exponential. (c) Expand its average in a series of cummulants. (d) Go back to the original representation. The final result is:

$$\xi \xi ^T(t)=e^{t\widehat{𝚲}}\xi _0\xi _0^T,$$
(12)
where $`\widehat{𝚲}`$ is a time-independent superoperator given by the perturbative expansion:

$`\widehat{𝚲}\widehat{𝐀}_0+\widehat{𝐀}_1`$ (13)
$`+{\displaystyle _0^{\mathrm{}}}d\tau \delta \widehat{𝐀}_1(t)e^{\tau \widehat{𝐀}_0}\delta \widehat{𝐀}_1(t\tau )e^{\tau \widehat{𝐀}_0}+\mathrm{},`$ (14)
with

$$\delta \widehat{𝐀}_1(t)=\widehat{𝐀}_1(t)\widehat{𝐀}_1.$$
(15)
Let $`L_{\mathrm{max}}`$ be the eigenvalue of $`\widehat{𝚲}`$ which has the largest real part. We find that the largest Lyapunov exponent $`\lambda `$ is related to the real part of $`L_{\mathrm{max}}`$:

$$\lambda =\frac{1}{2}\text{Re}\left(L_{\mathrm{max}}\right).$$
(16)
In Eq. (13) we give explicitly only the first two cummulants, the dots stand for third cummulants and higher order ones. The perturbative parameter can be understood as the product of two quantities. The first one, let’s call it $`\sigma `$, characterizes the amplitude of the fluctuations of $`\delta \widehat{𝐀}_1(t)`$. The second, $`\tau _c`$, is a typical (the largest relevant) correlation time of $`\delta \widehat{𝐀}_1(t)`$. Thus, the second cummulant is of the order of $`\sigma ^2\tau _c`$, the third one is of the order of $`\sigma ^3\tau _c^2`$, and so on. If all cummulants were summed up, Eq. (12) would be exact in the long-time regime $`t\tau _c`$ .

From now on, we restrict our analysis to the propagator $`\widehat{𝚲}`$ truncated at the second order, i.e., Eq. (13) without the dots. This approximation will be better the smaller $`\sigma \tau _c`$. However, if $`\widehat{𝐀}_1(t)`$ is not far from a Gaussian process, the validity of the second order approximation may extend outside the perturbative region $`\sigma \tau _c1`$. In the exceptional case that $`\widehat{𝐀}_1(t)`$ is a Gaussian process, cummulants higher than the second one will be strictly zero and the truncation will introduce no error.

To proceed further one needs the matrix of $`\widehat{𝚲}`$ in some basis. So, let us calculate $`\widehat{𝚲}𝐌`$, with $`𝐌`$ a symmetric matrix (it is easy to see that the truncation has not spoiled the symmetry of the density matrix). First notice that the exponentials of $`\widehat{𝐀}_0`$ represent no problem as they are finite polinomials,

$$e^{\tau \widehat{𝐀}_0}𝐐=[𝟙+\tau 𝐀_\mathrm{𝟘}]𝐐[𝟙+\tau 𝐀_\mathrm{𝟘}^𝕋],$$
(17)
for any matrix $`𝐐`$. Inserting this expression into (14) we arrive at

$`\widehat{𝚲}𝐌=(𝐀_0+𝐀_1)𝐌`$ (18)
$`+{\displaystyle _0^{\mathrm{}}}d\tau \delta 𝐀_1(t)[\delta \stackrel{~}{𝐀}_1(t\tau )𝐌+𝐌\delta \stackrel{~}{𝐀}_1^T(t\tau )]`$ (19)
$`+(\mathrm{})^T,`$ (20)
where $`(\mathrm{})^T`$ means “the previous terms transposed”, and

$$\delta \stackrel{~}{𝐀}_1(t\tau )=[𝟙+\tau 𝐀_\mathrm{𝟘}]\delta 𝐀_\mathrm{𝟙}(𝕥\tau )[𝟙\tau 𝐀_\mathrm{𝟘}].$$
(21)
Substituting (11) into (20) we arrive at the final result of the second-order perturbative approach

$`\widehat{𝚲}𝐌=\left(\begin{array}{cc}0& \\ 𝐕& 0\end{array}\right)𝐌`$ (22)
$`+{\displaystyle _0^{\mathrm{}}}𝑑\tau \left(\begin{array}{cc}0& 0\\ \tau & \tau ^2\end{array}\right)\left(\begin{array}{cc}\delta 𝐕\delta 𝐕^{}& 0\\ 0& \delta 𝐕\delta 𝐕^{}\end{array}\right)𝐌`$ (23)
$`+{\displaystyle _0^{\mathrm{}}}𝑑\tau \left(\begin{array}{cc}0& 0\\ \delta 𝐕& 0\end{array}\right)𝐌\left(\begin{array}{cc}\delta 𝐕^{}& 0\\ 0& \delta 𝐕^{}\end{array}\right)\left(\begin{array}{cc}\tau & \\ \tau ^2& \tau \end{array}\right)`$ (24)
$`+(\mathrm{})^T.`$ (25)
To abbreviate the notation, we have written $`\tau ^n`$ instead of $`\tau ^n𝟙`$; $`\delta 𝐕`$ and $`\delta 𝐕^{}`$ substitute $`\delta 𝐕(t)`$ and $`\delta 𝐕(t\tau )`$, respectively.

The largest Lyapunov exponent is buried into Eq. (22). To get an explicit expression one must diagonalize the matrix of $`\widehat{𝚲}`$. The outcome will be $`\lambda `$ as a function of the first two cummulants of the stochastic process $`𝐕(t)`$, i.e., averages and (integrated) two-times correlation functions:

$$V_{ij};_0^{\mathrm{}}𝑑\tau \tau ^n\delta V_{ij}(0)\delta V_{kl}(\tau ),n=0,1,2.$$
(26)
At first sight it may be thought that, as $`\widehat{𝚲}`$ is a superoperator, the matrix one should diagonalize is of the order of $`N^2\times N^2`$, then straightforward diagonalization would be out of the question for large $`N`$. Notwithstanding, $`\widehat{𝚲}`$ is an averaged object, and, as such, it possesses some symmetries which can be exploited to reduce the dimensionality of the problem to tractable levels, say $`N\times N`$. This will be illustrated with an example in Sect. IV. So, if desired, the largest Lyapunov exponent could be found by numerical diagonalization, at least for systems with $`N1000`$ degrees of freedom (provided one can estimate the correlation functions).

An alternative to exact diagonalization is approximate diagonalization, i.e., the diagonalization of the restriction of $`\widehat{𝚲}`$ to some small-dimension subspace. Notice that the problem we are dealing with is not very different from finding the ground state energy of a quantum system. However, in the quantum problem, the operator one must diagonalize is Hermitian, and it is well known that diagonalization in a truncated basis produces an upper bound to the exact ground-state energy. It frequently happens that this bound is close to the exact result, even if the ground-state wavefunction is not. In spite of the operator $`\widehat{𝚲}`$ not being Hermitian (see Sect. III) we still expect that diagonalization in a small basis will give a lower bound for the Lyapunov exponent. If the basis is suitably chosen, this estimate may be close to the result of the exact diagonalization.

To proceed with the construction of a basis for $`\widehat{𝚲}`$, we take advantage of the fact that $`\lambda `$ is independent of $`\xi _0`$, and simplify Eq. (12) further by averaging over an orthonormal set of initial tangent vectors, obtaining

$$\xi \xi ^T(t)=\frac{1}{2N}e^{t\widehat{𝚲}}𝟙.$$
(27)
This second averaging allows us to consider, instead of $`\widehat{𝚲}`$ itself, the restriction of $`\widehat{𝚲}`$ to the subspace spanned by the matrices $`\widehat{𝚲}^k𝟙`$, $`k=0,1,2,\mathrm{}`$. A look at the first terms of these sequence gives a hint for constructing an appropriate basis. The first term is the identity, the second one is

$`\widehat{𝚲}𝟙`$ $`=`$ $`\left(\begin{array}{cc}0& 𝟙𝐕\\ 𝟙𝐕& 0\end{array}\right)`$ (29)
$`+2{\displaystyle _0^{\mathrm{}}}𝑑\tau \left(\begin{array}{cc}0& \tau \delta 𝐕\delta 𝐕^{}\\ \tau \delta 𝐕\delta 𝐕^{}& (1\tau ^2)\delta 𝐕\delta 𝐕^{}\end{array}\right),`$
and so on.

## III The isotropic approximation

Typically, the diagonal elements of $`𝐕(t)`$ will be larger than the off-diagonal ones. This is evident in the case of translational invariance, where one has the property

$$V_{ii}=\underset{ji}{}V_{ij}.$$
(30)
Introducing a matrix $`𝐘`$ having all entries equal to one, i.e.,

$$𝐘_{ij}=1,i,j.$$
(31)
we can rewrite Eq. (30) as

$$\mathrm{𝐘𝐕}=\mathrm{𝐕𝐘}=0.$$
(32)
Then it is clear that Eq. (30) is also satisfied by $`𝐕`$, $`\delta 𝐕\delta 𝐕^{}`$, and by the higher moments of $`𝐕`$ that will appear in the blocks of $`\widehat{𝚲}^k𝟙`$, for $`k>1`$. So, in a first (crude) approximation one may be tempted to discard the off-diagonal elements of the moments of $`𝐕`$. If we also assume that all coordinates $`q_i`$ are statistically equivalent, and remind that the matrices $`\widehat{𝚲}^k𝟙`$ are symmetric, we arrive at the simplest approximation for diagonalizing $`\widehat{𝚲}`$. We call this approximation “isotropic”, and it consists in restricting $`\widehat{𝚲}`$ to the subspace spanned by the following three matrices:

$$𝐈_1=\left(\begin{array}{cc}& 0\\ 0& 0\end{array}\right),𝐈_2=\left(\begin{array}{cc}0& 0\\ 0& \end{array}\right),𝐈_3=\left(\begin{array}{cc}0& \\ & 0\end{array}\right).$$
(33)
These matrices are mutually orthogonal with respect to the standard Euclidean scalar product, i.e.,

$$\mathrm{Tr}(𝐈_i𝐈_j^T)\delta _{ij}.$$
(34)
Then the matrix elements of $`\widehat{𝚲}`$ with respect to the basis $`\{𝐈_1,𝐈_2,𝐈_3\}`$ are

$$\mathrm{\Lambda }_{ij}^{II}=\frac{\mathrm{Tr}([\widehat{𝚲}𝐈_j]𝐈_i^T)}{\mathrm{Tr}(𝐈_i𝐈_i^T)}.$$
(35)
Using Eq. (22) and skipping some simple algebra, we arrive at the 3$`\times `$3 matrix

$$𝚲^{II}=\left(\begin{array}{ccc}0& 0& 2\\ 2\sigma ^2\tau _c^{(1)}& 2\sigma ^2\tau _c^{(3)}& 2\mu \\ \mu +2\sigma ^2\tau _c^{(2)}& 1& 2\sigma ^2\tau _c^{(3)}\end{array}\right),$$
(36)
with the definitions

$`\mu `$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}𝐕,`$ (37)
$`\sigma ^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}\left(\delta 𝐕\right)^2,`$ (38)
$`\tau _c^{(k+1)}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\tau \tau ^kf(\tau ),`$ (39)
where we have introduced the normalized correlation function $`f(\tau )`$

$`f(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{N\sigma ^2}}\mathrm{Tr}\delta 𝐕(0)\delta 𝐕(\tau )`$ (40)
$`=`$ $`{\displaystyle \frac{1}{N\sigma ^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\delta V_{ij}(0)\delta V_{ij}(\tau ).`$ (41)
It is evident from Eq. (36) that the operator $`\widehat{𝚲}`$ is not Hermitian. Normalization of the basis $`𝐈_j`$ will not make $`𝚲^{II}`$ symmetric.

In the isotropic approximation the Lyapunov exponent is expressed in terms of the set of four parameters $`\mu `$ and $`\sigma ^2\tau _c^{(k+1)}`$, $`k=0,1,2`$. The parameters $`\mu `$ and $`\sigma `$ are, respectively, the mean and variance of the stochastic process $`𝐕(t)`$, and can in principle be obtained analytically by calculating the corresponding microcanonical averages. (In practice the calculations can be done in the canonical ensemble, and then connected with the microcanonical results by the formula of Lebowitz, Percus and Verlet .)

The characteristic time $`\tau _c^{(1)}`$ is naturally interpreted as the correlation time of the process $`𝐕(t)`$. Its calculation requires the knowledge of the autocorrelation functions of $`V_{ij}(t)`$, which are system dependent. Moreover, the correlation functions of a given system will in general depend on energy. However, if the functional form of the correlation function $`f(\tau )`$ is known (or conjectured), its parameters can also be calculated as thermal averages. For instance, if $`f(\tau )`$ is approximately Gaussian,

$$f(\tau )e^{\gamma \tau ^2},$$
(42)
the expansion of $`𝐕(0)𝐕(\tau )`$ around $`\tau =0`$ gives an explicit formula for the correlation time, namely

$$\frac{1}{\tau _c^{(1)}}=\left[\frac{2}{\pi \sigma ^2N}\text{Tr}\left(\frac{d𝐕}{dt}\right)^2\right]^{1/2}.$$
(43)
In this case $`\tau _c^{(2)}`$ and $`\tau _c^{(3)}`$ are trivially related to $`\tau _c^{(1)}`$:

$`\tau _c^{(2)}`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}\left[\tau _c^{(1)}\right]^2,`$ (44)
$`\tau _c^{(3)}`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}\left[\tau _c^{(1)}\right]^3.`$ (45)
A purely numerical calculation of $`\tau _c^{(k)}`$ may be very difficult, if not impossible, because correlation functions estimated from finite-length time series usually fail to damp as expected . Perhaps, a more sensible approach to the estimation of $`\tau _c^{(k)}`$ should start with a numerical study of the correlation functions; then a functional form for $`f(\tau )`$ could be proposed, based on the short-time behavior of the numerical correlation functions; finally, the parameters defining $`f(\tau )`$ would be calculated as suitable thermal averages. An alternative, more powerful approach, involves the use of “memory functions” : They are related in a one-to-one way to correlation functions and seem to be more amenable to simple approximations (see, e.g., Ref. ).

## IV Examples

In this section we analize the application of the perturbative theory of Sect. II to some simple models. We remark that, in principle, the theory is expected to be successful only in regimes where the Lyapunov exponent is very small. All the systems considered below exhibit regimes with vanishingly small Lyapunov exponents. It is understood that our discussion will be restricted to such regimes.

### A Mean field XY-Hamiltonian

Let us begin by analizing one special case in which the isotropic approximation of Sect. III is exact. Consider the one-dimensional Hamiltonian

$$H_1=\frac{1}{2}\underset{i=1}{\overset{N}{}}L_i^2+\frac{1}{2N}\underset{i,j=1}{\overset{N}{}}[1\mathrm{cos}(\theta _i\theta _j)].$$
(46)
This is the so-called mean-field $`XY`$-Hamiltonian. It represents a lattice of classical spins with infinite-range interactions. Each rotator is restricted to the unit circle and it is therefore described by an angle $`0<\theta _i2\pi `$ and its conjugate angular momentum $`L_i`$, with $`i=1,\mathrm{},N`$. At the critical energy $`E_c=3N/4`$ there is a second-order phase transition separating a disordered regime ($`E>E_c`$) from an ordered one ($`E<E_c`$).

In both limits $`E0,\mathrm{}`$ the Lyapunov exponent goes to zero. For a fixed energy $`E>E_c`$, $`\lambda `$ also goes to zero when $`N\mathrm{}`$. This behaviour has also been observed in a metastable disordered phase with $`E<E_c`$. The perturbative approach should be a good approximation in these regimes. Moreover, we argue that the infinite-range interactions justify the isotropic approximation.

All single-particle averages are equal, and, given that the forces are independent of the distances between spins, all two-particle averages must also be equal. So, one has

$`V_{ii}`$ $`=`$ $`c_1,i,`$ (47)
$`V_{ij}`$ $`=`$ $`c_2,ij.`$ (48)
Notice that translational symmetry, Eq. (30), implies that

$$c_2=\frac{c_1}{N1}.$$
(49)
This is the reason why the isotropic approximation will work in this case, i.e., off-diagonal matrix elements are indeed smaller than diagonal ones . But let us keep the discussion quantitave, and rewrite Eq. (47) as

$$𝐕=c_1𝟙+𝕔_\mathrm{𝟚}(𝐘𝟙),$$
(50)
with $`𝐘`$ defined in Eq. (31). Using the time-reversal symmetry of the stochastic process $`𝐕(t)`$, one can also show that

$$\delta 𝐕(0)\delta 𝐕(\tau )=c_1^{}𝟙+𝕔_\mathrm{𝟚}^{}(𝐘𝟙).$$
(51)
Then all blocks of $`\widehat{𝚲}𝟙`$ \[Eq. (29)\] belong to the subspace spanned by $`𝟙`$ and $`𝐘`$. Taking into account Eq. (32) and

$$𝐘^2=N𝐘,$$
(52)
we conclude that the blocks of all the sequence $`\widehat{𝚲}^k𝟙`$ belong to the subspace $`\{𝟙,𝐘\}`$. Thus, the relevant subspace for diagonalizing $`\widehat{𝚲}`$ is six-dimensional. It is spanned by $`𝐈_1,𝐈_2,𝐈_3`$ \[Eq. (33)\] and $`𝐘_1,𝐘_2,𝐘_3`$, with the definitions

$$𝐘_1=\left(\begin{array}{cc}𝐘& 0\\ 0& 0\end{array}\right),𝐘_2=\left(\begin{array}{cc}0& 0\\ 0& 𝐘\end{array}\right),𝐘_3=\left(\begin{array}{cc}0& 𝐘\\ 𝐘& 0\end{array}\right).$$
(53)
However, it can be shown (see Appendix) that the largest eigenvalue of the corresponding 6$`\times `$6 matrix coincides with that of the isotropic 3$`\times `$3 matrix, up to corrections of order $`1/N`$. In this way we have proven the validity of the isotropic 3$`\times `$3 approximation for one-dimensional systems with infinite-range forces.

### B Dilute gases

Consider now a one-dimensional gas with Hamiltonian

$$H_2=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2+\underset{i,j=1}{\overset{N}{}}\nu (q_iq_j),$$
(54)
where $`q`$ and $`p`$ are linear coordinates, and we assume that the potential $`\nu `$ is bounded. For large enough energies, this system is disordered and weakly chaotic. All particles, and all pairs of particles, are statistically equivalent. Then this problem is formally equivalent to the infinite range $`XY`$-Hamiltonian: The isotropic 3$`\times `$3 approximation becomes exact.

Of course, the statistical equivalence also holds for a dilute three-dimensional gas with short-range interactions. In this case, Barnett et al. have shown that the largest Lyapunov exponent is found by diagonalizing a 4$`\times `$4 matrix .

### C $`\alpha XY`$-Hamiltonian

There are cases in which no strong reasons exist to believe that the isotropic approximation will work satisfactorily. Consider, for instance, the arbitrary-range analog of the $`XY`$-Hamiltonian :

$$H_3=\frac{1}{2}\underset{i=1}{\overset{N}{}}L_i^2+\frac{1}{2\stackrel{~}{N}}\underset{i,j=1(ij)}{\overset{N}{}}\frac{1\mathrm{cos}(\theta _i\theta _j)}{r_{ij}^\alpha }.$$
(55)
The parameter $`\alpha `$ sets the range of the interactions: $`\alpha =0`$ recovers the mean-field case, $`\alpha =\mathrm{}`$ corresponds to first-neighbors couplings. The prefactor $`\stackrel{~}{N}`$ (a function of $`N`$ and $`\alpha `$) is included to make the system “pseudo-extensive” . Periodic boundary conditions are assumed, and $`r_{ij}`$ is the minimum between $`|ij|`$ and $`N|ij|`$. For any value of $`\alpha `$ there exist (i) a low-energy regime of harmonic oscillators weakly coupled by non-linear forces and (ii) a high-energy disordered phase where the spins rotate almost freely. We expect our theory to produce good estimates for the largest Lyapunov exponent in both low- and high-energy regimes.

If the forces are not infinite-range but just long-ranged, the isotropic approximation will still give good estimates in weakly chaotic regimes. Evidence supporting this statement can be found in Refs.  (geometric method) and  (random matrix approach), where some kind of “isotropic” approximations were used to derive scaling laws for $`\lambda `$ in the high-energy regime, in good agreement with numerical simulations .

For $`\alpha `$ not too small, it may be necessary to improve the isotropic approximation by diagonalizing $`\widehat{𝚲}`$ in a larger basis. In this case, the statistical equivalence holds for all pairs of particles separated by the same distance. This means that the blocks of $`\widehat{𝚲}^k𝟙`$ are symmetric and cyclical, i.e., the matrix elements only depend on the distance $`r_{ij}`$. A basis can be constructed starting from the $`N`$$`\times `$$`N`$ matrix $`𝐒`$ of a cyclical shift:

$$S_{ij}=\delta _{i,j+1},$$
(56)
where it is understood that $`j+1`$ must be taken modulo $`N`$. Then the set of symmetrical matrices

$$𝚺_k𝐒^k+𝐒^k,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}kN/2,$$
(57)
is a basis for the blocks of $`\widehat{𝚲}^m𝟙`$. A suitable basis for diagonalizing $`\widehat{𝚲}`$ is the set

$$\left(\begin{array}{cc}𝚺_i& 0\\ 0& 0\end{array}\right),\left(\begin{array}{cc}0& 0\\ 0& 𝚺_j\end{array}\right),\left(\begin{array}{cc}0& 𝚺_k\\ 𝚺_k& 0\end{array}\right),$$
(58)
with $`0i,j,kN/2`$. The length of this basis is 3$`N`$/2. Notwithstanding, we expect that a small subset of this basis will be enough to get a satisfactory convergence to the largest eigenvalue of $`\widehat{𝚲}`$. Even in the worst case of no truncation at all, numerical diagonalization is possible for relatively large systems.

## V Connection with the geometric method

In Sect. III we motivated the isotropic approximation by arguing that, in a first approach, one can neglect the off-diagonal matrix elements of the blocks of $`\widehat{𝚲}^k𝟙`$. Then, in Sect. IV, we proved that this approximation is indeed justified in various cases. Looking back to the results of Sects. II and III, we realize that the isotropic approximation is equivalent to postulating an “effective” system of equations,

$$\dot{\xi }_i=\left(\begin{array}{cc}0& 1\\ K(t)& 0\end{array}\right)\xi _i,$$
(59)
where $`\xi _i=(\delta q_i,\delta p_i)`$ is the projection of the tangent vector $`\xi `$ on the subspace of the $`i`$-th degree of freedom. The equations above represent the evolution of a typical component of $`\xi `$. In this sense they could alternatively be called “mean-field” or “single-particle” equations. The scalar object $`K(t)`$ is an effective random process which substitutes the Hessian $`𝐕(t)`$, and is in principle unknown. However, its first two cummulants can be identified in the following way. First solve Eq. (59) for the average of $`\xi _i^2`$ by using second order perturbation theory, as done in Sect. II \[just change $`𝐕`$ by $`K`$, and set $`N`$=1 in Eq. (22)\]. Notice that, as $`K(t)`$ is a real number, the blocks of the effective $`\widehat{𝚲}`$ are also real numbers, and the “isotropic approximation” is exact now. Then the matrix one must diagonalize to obtain the Lyapunov exponent is exactly that of Eq. (36), provided one makes the identifications:

$`K`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}𝐕,`$ (60)
$`\delta K(0)\delta K(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}\delta 𝐕(0)\delta 𝐕(\tau ),`$ (61)
with $`\delta K=KK`$. From this point of view, the isotropic basis $`\{𝐈_1,𝐈_2,𝐈_3\}`$ is a single-particle basis. It is the most natural one in the sense that it treats all degrees of freedom on the same footing.

The perturbative-isotropic approximation, as presented above, is very similar to the geometric approach. In fact, in the geometric method  an effective equation like (59) is proposed, containing the unknown process $`K^{}(t)`$. Then, it is argued that the first two cummulants of $`K^{}`$ are related to the Laplacian of the potential, $`\mathrm{}𝒱(t)`$:

$`K^{}`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{}𝒱,`$ (62)
$`\delta K^{}(0)\delta K^{}(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left(\delta \mathrm{}𝒱\right)^2\overline{\tau }\delta (\tau ),`$ (63)
where $`\delta \mathrm{}𝒱=\mathrm{}𝒱\mathrm{}𝒱`$, and $`\overline{\tau }`$ is the correlation time of the process $`K^{}(t)`$, which is assumed to be delta-correlated.

It is obvious that the averages of both processes $`K`$ and $`K^{}`$ coincide because $`\mathrm{Tr}𝐕=\mathrm{}𝒱`$. So, the differences between the geometric method and the perturbative-isotropic approach appear only in the fluctuations. The similarity between both theories could be enhanced by relaxing the delta-correlation assumption of the geometric method, and substituting Eq. (63) by

$$\delta K^{}(0)\delta K^{}(\tau )=\frac{1}{N}\delta \mathrm{}𝒱(0)\delta \mathrm{}𝒱(\tau ).$$
(64)
But even so, we have not been able to find any analytical relationship between the correlation functions of $`K`$ and $`K^{}`$: In principle, the difference between both is non-negligible, and both effective theories will lead to different estimates for the Lyapunov exponent. We expect that numerical simulations will decide which estimate is better.

One comment about the correlation time $`\overline{\tau }`$ of Eq. (63) is in order: Geometric arguments lead to the estimate

$$\overline{\tau }=\frac{\pi \sqrt{\overline{\mu }}}{2\sqrt{\overline{\mu }(\overline{\mu }+\overline{\sigma })}+\pi \overline{\sigma }},$$
(65)
with $`\overline{\mu }K^{}`$ and $`\overline{\sigma }^2(\delta K^{})^2`$. Some slightly different expressions have also been proposed . The criterion for testing the accuracy of theses estimates has been the agreement between the geometric estimate for $`\lambda `$ and numerical simulations, i.e., the goodness of fit (which is indeed excellent in some cases). To our knowledge, there is no independent test of the expression (65), or others, in the literature. Accordingly, a precise definition of $`\overline{\tau }`$ seems to be lacking. (Is $`\overline{\tau }`$ equal to the integral of the normalized autocorrelation function of $`\mathrm{}𝒱(t)`$?) This is a point which affects the consistency of the geometric method. Unless a definition is given, to some extent, $`\overline{\tau }`$ will have the status of a fitting parameter. Comparisons of the geometric method with other theories will have to take this fact into account.

## VI Summary

We showed that the evolution equation in tangent space can be thought of as a stochastic differential equation with multiplicative noise. Then, an analytical estimate for the largest Lyapunov exponent of a many-particle system in equilibrium was derived by using standard perturbative techniques. Our analysis has been focused on the second-order approximation. In this case the Lyapunov exponent can be obtained by diagonalizing a matrix whose entries are calculated from the first two cummulants of the Hessian of the potential energy, i.e., the averages $`V_{ij}`$ and the correlation functions $`\delta V_{ij}(0)\delta V_{kl}(\tau )`$. The dimension of this matrix is, in principle, of the order of $`N`$$`\times `$$`N`$, but we have proposed the conjecture, based on an analogy with the Hermitian problem, that diagonalization in a truncated basis may be enough to obtain satisfactory results.

In the crudest approximation, which consists in choosing the three-dimensional isotropic basis of Eq. (33), the Lyapunov exponent is extracted from a 3$`\times `$3 matrix. We argued that this “isotropic approximation” is equivalent to modeling the tangent dynamics of the many-particle system by an “effective” process $`K(t)`$ for a single degree of freedom. In this way we established a connection with the so-called geometric method, the alternative effective theory for estimating the Lyapunov exponent. Both theories are very similar, but differ at the point of the definition of the correlation function of $`K(t)`$. The difference is non-trivial and is expected to lead to different predictions.

In special cases, e.g., one-dimensional lattice-systems with infinite-range interactions, we have been able to prove that the isotropic approximation is exact. However, in the general case, it may be necessary to consider larger bases. We have given examples where these bases are constructed by following the symmetries of the moments of $`𝐕(t)`$.

The theory we have presented is perturbative. Loosely speaking, we expect to obtain good estimates of Lyapunov exponents in weakly chaotic regimes. More quantitatively, the domain of validity of the theory is controlled by the “Kubo number” $`\sigma \tau _c`$, which quantifies the strength of the fluctuations $`\delta 𝐕(t)`$. For a given system, it is difficult to say a priori in which regimes the theory will work satisfactorily. This question and others, like the validity of the isotropic approximation, and its comparison with the geometric method, will be decided with the aid of forthcomming numerical simulations.

## Acknowledgements

We are grateful to A. O. Caldeira, C. H. Lewenkopf, F. Nobre, A. M. Ozorio de Almeida, H. M. Pastawski, A. Robledo, A. M. C. de Souza, and C. Tsallis for fruitful discussions. We acknowledge Brazilian Agencies CNPq, FAPERJ and PRONEX for financial support.

## A The infinite-range case

We have seen in Sect. IV that in the case of a one-dimensional system with infinite range interactions the subspace spanned by the matrices $`\widehat{\mathrm{\Lambda }}^k𝟙`$ is six-dimensional. An ortoghonal basis for this subspace is the set $`\{𝐈_1,𝐈_2,𝐈_3,𝐙_1,𝐙_2,𝐙_3\}`$, where

$$𝐙_k=𝐘_k𝐈_k.$$
(A1)
The 6$`\times `$6 matrix of $`\widehat{𝚲}`$ can be naturally split into four blocks of size 3$`\times `$3. The block $`I`$-$`I`$ has already been calculated \[Eq. 36\]. Let’s now calculate the block $`I`$-$`Z`$, i.e.,

$$\mathrm{\Lambda }_{ij}^{IZ}=\frac{\mathrm{Tr}(\widehat{𝚲}𝐙_j)𝐈_i}{\mathrm{Tr}𝐈_i^2}.$$
(A2)
By setting $`𝐕=0`$ in Eq. (22) we obtain the operator $`\widehat{𝚲}_0`$ ($`\widehat{𝚲}`$ in the absence of interactions). It has the following properties:

$`\widehat{𝚲}𝐘_j`$ $`=`$ $`\widehat{𝚲}_0𝐘_j,`$ (A3)
$`\mathrm{Tr}\left(\widehat{𝚲}_0𝐘_j\right)𝐈_i`$ $`=`$ $`\mathrm{Tr}\left(\widehat{𝚲}_0𝐈_j\right)𝐈_i.`$ (A4)
Using these two properties together with Eq. (A1) one arrives at

$$\mathrm{\Lambda }_{ij}^{IZ}=\mathrm{\Lambda }_{0,ij}^{II}\mathrm{\Lambda }_{ij}^{II}.$$
(A5)
Analogously one obtains

$`\mathrm{\Lambda }_{ij}^{ZI}`$ $``$ $`{\displaystyle \frac{1}{N}}\mathrm{\Lambda }_{ij}^{IZ},`$ (A6)
$`\mathrm{\Lambda }_{ij}^{ZZ}`$ $``$ $`\mathrm{\Lambda }_{0,ij}^{II}+{\displaystyle \frac{1}{N}}\mathrm{\Lambda }_{ij}^{II},`$ (A7)
where the symbol $``$ means that terms of relative size $`1/N`$ have been discarded. The 6$`\times `$6 matrix reads

$$𝚲\left(\begin{array}{cc}𝚲^{II}& 𝚲_1^{II}\\ \frac{1}{N}𝚲_1^{II}& 𝚲_0^{II}+\frac{1}{N}𝚲_1^{II}\end{array}\right),$$
(A8)
with the definition $`𝚲_0^{II}+𝚲_1^{II}=𝚲^{II}`$. Then it can be checked that the matrix above has three zero eigenvalues while the remaining three are the eigenvalues of the matrix

$$𝚲_0^{II}+\frac{N+1}{N}𝚲_1^{II}𝚲^{II}.$$
(A9)
Thus the isotropic approximation is essentially exact for the infinite-range $`XY`$-Hamiltonian.


format: text

text:
# X-Ray Spectra of Intermediate-Luminosity, Radio-Loud Quasars

## 1 Introduction

A major outstanding question in our ongoing efforts to understand the central engines of active galactic nuclei (AGNs) is the origin of the dichotomy between radio-loud and radio-quiet objects. Although both radio-loud and radio-quiet AGNs are thought to be powered by accretion onto a supermassive black hole, it is not known why AGNs in the former class can produce collimated, relativistic radio jets, which power the giant radio lobes at large distances from their host galaxies, while AGNs in the latter class cannot (see Livio, 1997, for a review). One possible explanation is that the central engines of the two classes of objects are the same, but systematic differences in the properties of the host galaxies are the cause of the observed dichotomy (e.g., Blandford & Levinson, 1995; Fabian & Rees, 1995). Another possibility is that there is a systematic difference between the central engines of radio-loud and radio-quiet AGNs, which is intimately related to their ability to form jets. Such a difference could be the size of the accretion disk (invoked in some wind models; e.g., Blandford & Payne, 1982), the spin of the black hole (thought to be the ultimate power source of the jets; e.g., Blandford & Znajek, 1977; Wilson & Colbert, 1995; Meier, 1999), or the structure of the inner accretion disk (an ion torus or advection-dominated accretion flow, a.k.a. ADAF, which allows some of the accreting matter to form an outflow; e.g., Rees et al., 1982; Narayan & Yi, 1995; Blandford & Begelman, 1999).

X-ray observations are a means by which one can investigate the last two of the above possibilities, since the X-rays are produced in the inner regions of the accretion flow. In particular, one can look for systematic differences in the X-ray spectral properties of radio-loud and radio-quiet AGNs that may indicate differences in the structure of the inner accretion disk or the spin of the central black hole. Thus we and other groups have been carrying out systematic studies of the X-ray properties of radio-loud AGNs using archival as well as proprietary ASCA and RXTE observations. These studies have confirmed previously suspected systematic differences between radio-loud and radio-quiet AGNs (Zdziarski et al., 1995; Woźniak et al., 1998; Sambruna, Eracleous, & Mushotzky, Sambruna et al.1999; Eracleous, Sambruna, & Mushotzky, 2000). In particular, two spectral features that are the hallmarks of the X-ray spectra of Seyfert galaxies (e.g., Nandra & Pounds, 1994; Nandra et al., 1997b): the fluorescent Fe K$`\alpha `$ line at 6.4 keV and the Compton “reflection” hump at energies greater than 10 keV, are considerably weaker in the spectra of broad-line radio galaxies (hereafter BLRGs; Zdziarski et al., 1995; Woźniak et al., 1998; , Sambruna et al.1999; Eracleous et al., 2000). The above observational results suggest that the structure of the inner accretion disk of BLRGs is different from that of Seyfert galaxies, although the nature of this difference remains to be understood in detail. Another potentially important trend is the X-ray Baldwin effect (Nandra et al., 1997c), according to which the Fe K$`\alpha `$ lines of luminous quasars are systematically weaker (i.e., of lower equivalent width; hereafter EW) than those of the lower-luminosity Seyfert galaxies. Because the sample of objects studied by Nandra et al. (1997c) includes a mix of radio-loud and radio-quiet quasars (hereafter QSRs and QSOs, respectively), it is not yet clear whether and how such an effect applies to radio-loud and radio-quiet AGNs separately. Finally, the immediate environment of the X-ray source appears to be different in radio-loud and radio quiet AGNs as evidenced by the properties of absorbing material along the line of sight. More specifically, the signature of a “warm” (i.e., ionized) absorber, seems to be ubiquitous in Seyfert galaxies (e.g., Reynolds, 1997; George et al., 1998) but rare in broad-line radio galaxies (, Sambruna et al.1999). Instead, radio-loud AGNs seem to display the signature of significant absorption by neutral matter in their X-ray spectra.

Here we present the results of new observations of three intermediate-luminosity QSRs (3C 323.1, PKS 2349–01, 4C 74.24; $`L_\mathrm{X}10^{45}\mathrm{erg}\mathrm{s}^1`$) with ASCA and BeppoSAX, which were carried out as part of our larger program. The targets were specifically selected to bridge the gap in luminosity between the nearby and well-studied BLRGs and the more luminous QSRs included in the sample of Nandra et al. (1997c). All three are associated with large, double-lobed radio sources (Antonucci, 1985; Kellerman et al., 1994; Riley et al., 1988), which ensures that our line of sight is not close to the axis of their radio jet, thus the X-rays we observe are not contaminated by beamed emission from the jet itself. The goal of the observations was to investigate the properties of intermediate-luminosity QSRs and to check whether the systematic differences between Seyfert galaxies and BLRGs extend to higher luminosity objects. In this respect, 4C 74.26 is a particularly interesting object because its ASCA X-ray spectrum (Brinkmann et al., 1998; , Sambruna et al.1999) showed features that are hallmarks of Seyfert galaxies but not BLRGs: a “warm” absorber with a column density of $`3.5\times 10^{21}\mathrm{cm}^2`$ and a significant Compton “reflection” hump. It was also the only QSR in the collection of (Sambruna et al.1999) with a detectable Fe K$`\alpha `$ line.

In §2 we describe the observations and basic data analysis while in §3 we present the results of fitting models to the spectra. In §4 we discuss our results and we present our final conclusions. Throughout this paper we adopt a Hubble constant of $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and a deceleration parameter of $`q_0=0.5`$.

## 2 Observations, Data Screening, and Timing Analysis

ASCA observations of PKS 2349–01 and 3C 323.1 were carried out in 1999 December and 2000 February, respectively. Table 1 reports the observation log. The SIS detectors were operated in 1-CCD FAINT mode half of the time and in 1-CCD BRIGHT mode for the remaining portion, while the GIS detectors were operated in PH mode. The data were screened within FTOOLS/XSELECT v.4.2 (Blackburn, Greene, & Pence, 1994; Ingham, 1994) following the procedures described in Eracleous & Halpern (1998) and Eracleous et al. (1996). The final net exposure times are included in Table 1. After screening, we extracted spectra and light curves for each of the $`\mathrm{𝐴𝑆𝐶𝐴}`$ detectors. The `FAINT` mode data were converted into `BRIGHT2` mode and combined to the `BRIGHT` mode data, to take advantage of the full exposure time and maximize the signal-to-noise ratio ($`S/N`$).

4C 74.26 was observed in 1999 May (Table 1) by the BeppoSAX narrow field instruments: the Low- and Medium-Energy Concentrator Spectrometers, and the Phoswich Detector System (hereafter LECS, MECS, and PDS, respectively). The data from two of the three MECS units were merged to increase the $`S/N`$. The exposure times resulting from the application of standard screening criteria to account for passage through the South-Atlantic Anomaly, Earth occultation, etc. (Fiore et al., 1999), are listed in Table 1. The source was detected in the PDS band up to energies of order 100 keV, with a count rate of 0.52 $`\pm `$ 0.04 c s-1 (Table 1). No bright sources were found in the LECS and MECS fields of view, thus enhancing our confidence that the source detected by the PDS is indeed the target quasar.

Figure 1 shows the ASCA and BeppoSAX light curves of the three targets. For PKS 2349–01 and 3C 323.1, we show the 2–10 keV light curves from the two ASCA GIS detectors combined, binned at 30-minute intervals. For 4C 74.26, we show the 2–10 keV light curve from the MECS instrument, binned at 30-minute intervals, as well as the 15–125 keV light curve from the PDS instrument, binned at 2-hour intervals. No systematic variability is evident in any of the light curves save for a brief, 25% excursion in the case of PKS 2349–01. To search for short-time scale fluctuations we applied the excess variance test (Nandra et al., 1997a; Turner et al., 1999) to the light curves binned at both 256-second and 1800-second intervals. The excess variance measured from the 256-second light curves is $`(5.4\pm 0.1)\times 10^2\mathrm{s}^1`$, $`(5.5\pm 0.3)\times 10^2\mathrm{s}^1`$, and $`(4.2\pm 0.1)\times 10^2\mathrm{s}^1`$ for PKS 2349–01, 3C 323.1, and 4C 74.26, respectively. These values are higher than those of Seyfert galaxies of comparable luminosity (see, for example Leighly, 1999; Turner et al., 1999) but they are similar to the values found for narrow-line Seyfert 1 galaxies of comparable luminosity. The corresponding values of the excess variance for the 1800-second light curves are $`(5.8\pm 0.1)\times 10^2\mathrm{s}^1`$, $`(3.7\pm 0.1)\times 10^2\mathrm{s}^1`$, and $`(1.8\pm 0.1)\times 10^2\mathrm{s}^1`$ for PKS 2349–01, 3C 323.1, and 4C 74.26, respectively. The observed fluctuations are unlikely to be related to emission from the jets of these objects because the radio morphology suggests that the the jets are oriented at large angles to the line of sight (Antonucci, 1985; Kellerman et al., 1994; Riley et al., 1988). This conclusion is also supported by the high-energy light curve of 4C 74.26, which is rather steady over the course of the observation, in marked contrast to the light curves of blazars.

## 3 Model Fits to the Spectra

Models were fitted to the ASCA and BeppoSAX spectra using XSPEC v.11.0.1 (Arnaud, 1996), in the energy ranges where the calibration is well established and the background is low (see Table 1). The spectra were rebinned such that each new bin contained at least 20 counts in order to validate the use of the $`\chi ^2`$ statistic. We used the most current versions of the response matrices and effective area curves (September 1, 1997 files for the BeppoSAX instruments; February 3, 1995 response matrices for the ASCA GIS, and all other ASCA calibration files generated specifically for our observations).

To check the self-consistency of the spectra obtained by different ASCA detectors, we first fitted the individual SIS and GIS spectra. We then fitted the same model to the spectra from all detectors simultaneously. To take account of any uncalibrated differences between the effective areas of the different instruments, we allowed the normalization of the corresponding spectra to vary freely. An analogous procedure was adopted for the spectra from the BeppoSAX detectors, as we discuss further below. As with the ASCA spectra, the normalization constants of spectra from different instruments were also allowed to vary, but the PDS normalization was constrained to be between 0.77 and 0.83 of the MECS normalization. The best-fitting value for the PDS/MECS ratio was found to be $`0.83_{0.06}^{+0}`$. The best-fitting LECS/MECS ratio, which was left free to vary, was $`0.70\pm `$0.02.

The goals of the spectral analysis are to determine the model that best describes the continuum, including possible excess absorption at low energies and the Compton hump at $`E\stackrel{>}{}10`$ keV, and check for the presence of the Fe K$`\alpha `$ emission line. To determine the best-fitting continuum models, we adopted the following procedure. At first, the spectra were fitted by excluding the rest energy range 5–7 keV, where the Fe K$`\alpha `$ line contributes most. Once the best-fitting continuum model was found, the 5–7 keV range was added back to the spectra and a Gaussian line was added to the model. This procedure ensures the continuum and the emission line remain unaffected by one another in the fitting process, an especially important concern if the line is broad (e.g., Nandra et al., 1997b).

While we investigated the presence of excess absorption at low energies in the three targets, a detailed analysis was possible only in the case of 4C 74.26, as the LECS sensitivity extends to lower energies than the SIS and the spectral calibration is more certain. For each ASCA detector we only use data in the energy range given in Table 1 in order to bypass calibration problems associated with the SIS efficiency loss at low energies (Yaqoob et al., 2000). Excluding the lowest-energy bins in the SIS spectra results in larger uncertainties in the determination of the absorbing column, however the column densities turn out to be high enough in both ASCA targets that they can be determined fairly well even by the GIS detectors. The presence of a Compton “reflection” hump can be determined with confidence only in the BeppoSAX data because of the higher-energy sensitivity of the PDS compared to the SIS and GIS detectors.

### 3.1 The Continuum

The best-fitting model for the ASCA spectra of PKS 2349–01 and 3C 323.1 is an absorbed power law with column density in excess of the Galactic value plus a Gaussian emission line. The best-fitting models, parameters, and corresponding error bars (at 90% confidence for two parameters of interest) are summarized in Table 2, while the spectra are shown in Figure 2 with models superposed. All the models include an absorption component, fixed to the Galactic column density, and a second absorption component (due to neutral gas) at the source’s redshift with the column density left free to vary in the fit. The implicit assumption of this procedure is that any excess absorption is at the source’s rest-frame. The simple power-law model also provides a good description of the LECS+MECS ($`E<10`$ keV) spectrum of 4C 74.26 but it leaves large residuals at high energies, in the portion of the spectrum sampled by the PDS, as shown in Figure 2c.

It is noteworthy that the absorbing column detected in the spectrum of 4C 74.26 is rather high, namely $`3\times 10^{21}\mathrm{cm}^2`$, a few times higher than what is measured for the other two objects. Analysis of the ASCA spectrum of this object by Brinkmann et al. (1998) suggested the presence of an intrinsic, ionized absorber (“warm” absorber). That result, however, was not confirmed in the independent analysis of the same data by (Sambruna et al.1999). To investigate this question further we fitted the LECS spectrum with two different models representing an ionized absorber: (a) a single power law with Galactic and intrinsic absorption, plus a redshifted absorption edge, and (b) an ionized absorber model according to Done et al. (1992, implemented by the XSPEC routine absori). In neither case did we find a significant improvement of the fit: the total value of $`\chi ^2`$ decreases because of the additional free parameters of these more sophisticated absorption models, but the improvement is significant only at the 50–65% level (according to the F-test). Our conclusion, therefore, is that although a warm absorber provides an acceptable fit, it is not required by the data.

The high-energy residuals in the simple power-law fit to the spectrum 4C 74.26 indicate the presence of an additional continuum component emerging at energies $`\stackrel{>}{}`$10 keV. The following two possibilities were considered in detail:

1. The high-energy excess flux is due to Compton “reflection” of the primary continuum by neutral, dense matter in the immediate vicinity of the X-ray source (Lightman & White, 1988; Guilbert & Rees, 1988). Thus, we fitted the BeppoSAX spectrum with a Compton “reflection” model (pexrav within XSPEC Magdziarz & Zdziarski, 1995). This model has several additional free parameters, which include the inclination angle, $`i`$, of the reprocessing slab relative to the line of sight, the solid angle, $`\mathrm{\Omega }`$, that the slab subtends to the primary X-ray source, and the upper cut-off energy of the primary power-law spectrum. In the particular case of 4C 74.26 we can use the radio properties to constrain the inclination angle of the radio jet (hence the inclination angle of the disk). The brightness contrast between the small-scale jet and counter-jet (observed by VLBI), if ascribed to relativistic beaming, yields an upper limit to the jet inclination angle of $`i49^{}`$ (Pearson et al., 1992). Requiring that the size of the large-scale radio source does not exceed the size of giant radio sources observed in radio galaxies can yield a lower limit on the inclination. Such giant radio source sizes typically range up to 2.0–2.5 Mpc (e.g., Palma et al., 2000; Lara et al., 2001). Assuming that the radio source in 4C 74.26 is not exceptionally large and adopting 2 Mpc as an upper limit to its true size, we obtain $`i37^{}`$. Thus we can express the inclination angle as $`i=43\pm 6`$ degrees. 111The constraints on the inclination angle of the jet derive from ratios of fluxes or lengths; thus they do not depend on $`H_0`$. In Table 2 we give the best fitting model parameters for the case where the inclination angle of the axis of the reprocessing slab was held fixed to the value inferred from the radio morphology. In Figure 4 we show how the inferred reprocessor solid angle depends on the value of the inclination by plotting the 68, 90, and 99% confidence contours in the $`\mathrm{cos}i\mathrm{\Omega }/2\pi `$ plane. For inclination angles between $`37^{}`$ and $`49^{}`$, which represent the range inferred from the radio properties, the solid angle is consistent with $`2\pi `$. The fit with the Compton “reflection” model is an improvement over the single power law at the 98.7% confidence level, according to the F-test.
2. The excess flux could indicate a contribution from the base of the jet. This was modeled by adding a second power law to the continuum, which did not yield a significantly different description, however, leading us to disfavor this option. More specifically, the spectral indices of the two power-law components converge to values that are indistinguishable from each other and the total $`\chi ^2`$ value is higher than the Compton “reflection” model, which has more free parameters. Moreover, emission from a jet is disfavored by the radio properties of this object (mainly the morphology), which suggest that the jet is oriented at a large angle to the line of sight ($`i>37^{}`$).

Although the ASCA bandpass is not well-suited to the detection of the Compton “reflection” continuum, we nevertheless searched for it in the spectra of PKS 2349–01 and 3C 323.1. The main reason for embarking on this exercise is that the adopted continuum model could have an influence on the measured properties of the Fe K$`\alpha `$ lines, which we study in the next section (see the discussion in the context of 4C 74.26 in §3.3). In the case of 3C 323.1 we find no evidence for a Compton “reflection” continuum. In fact, a Compton “reflection” model converges to a model with negligible “reflection,” which is in effect the original simple power-law model. In the case of PKS 2349–01 the Compton “reflection” model results in a slightly lower $`\chi ^2`$ value than the simple power-law model, but this apparent improvement is not significant (the F-test gives a 66% probability of improvement by chance). But the two continuum models yield significantly different properties for the Fe K$`\alpha `$ line, as we describe below. This effect has been known for some time, and has been discussed by a number of authors, including most recently Zdziarski & Grandi (2001).

### 3.2 The Fe K$`\alpha `$ Emission Line

To search for the presence of an Fe K$`\alpha `$ emission line, which is a ubiquitous feature of the X-ray spectra of Seyfert galaxies, we modeled its profile with a Gaussian. We fitted the continuum as described above and fixed the continuum parameters before adding the Gaussian line profile to the model. To evaluate whether the line is unambiguously detected and resolved we searched the 2-dimensional parameter space defined by the line photon flux and the energy dispersion of the Gaussian profile. The results of this search are shown in Figure 5 in the form of confidence contours in the line EW $`vs`$ FWHM plane (the EW and FWHM are directly proportional to the line photon flux and Gaussian energy dispersion, respectively). These results are also summarized in Table 3.

In the case of 4C 74.26 the line is clearly detected with a rest EW of $`170_{100}^{+180}`$ eV. The energy of the Fe K$`\alpha `$ emission line encompasses the range 6.4 to 6.9 keV within 90% confidence. Thus, we confirm the detection of the line in 4C 74.26 with ASCA, as reported by Brinkmann et al. (1998) and (Sambruna et al.1999). In the case of 3C 323.1 the line is only marginally detected, i.e., its EW is non-zero at 90% confidence but consistent with zero at 99% confidence (for an assumed energy of 6.4 keV). The best-fitting FWHM is about $`5\times 10^4\mathrm{km}\mathrm{s}^1`$, which is comparable with what is found in Seyfert galaxies by Nandra et al. (1997b).

In the case of PKS 2349–01, the measured properties of the Fe K$`\alpha `$ line depend sensitively on the continuum model. If we adopt a simple power-law model, the line is unambiguously detected at the 99% confidence level with a rest EW of $`230\pm 120`$ eV and a FWHM of $`6_3^{+6}\times 10^4\mathrm{km}\mathrm{s}^1`$. If, however, we adopt a model consisting of a power-law plus Compton “reflection,” fixing the solid angle of the reprocessor to $`\mathrm{\Omega }=2\pi `$, the line is about a factor of 2 weaker (see Table 3). If we fix the solid angle of the reprocessor to $`\mathrm{\Omega }=4\pi `$, its maximum physically reasonable value, the line is only marginally detected (at the 68% confidence level). The upper limit to its rest EW is 200 eV, assuming a FWHM of $`5\times 10^4\mathrm{km}\mathrm{s}^1`$. This ambiguity cannot be resolved with the data at hand. In order to measure the properties of the Fe K$`\alpha `$ lines with confidence, we need either a 2–10 keV spectrum of much higher $`S/N`$ or a spectrum of comparable $`S/N`$ extending to at least 30 keV, which can yield meaningful constraints on the properties of the Compton “reflection” component. Finally, we also note that the observed EWs of the Fe K$`\alpha `$ lines and the strength of the Compton “reflection” components are consistent with each other, albeit within the large uncertainties, as one would expect in the context of reprocessing models such as those of George & Fabian (1991). Such models attribute the origin of both the Fe K$`\alpha `$ line and the Compton “reflection” component to the same reprocessing medium and require that their strengths are related.

### 3.3 Comparison With Previous Observations of 4C 74.26 by ASCA

Our preferred model for 4C 74.26 consists of a power-law continuum modified by Compton “reflection” at high energies and by absorption at low energies. The line energy obtained from the model could not be tightly contrained, falling between 6.4 and 6.9 keV at a 90% confidence level. Within uncertainties, our results are largely consistent with the results of Brinkmann et al. (1998), who considered our preferred model among others. However, there is a systematic discrepancy between the properties of the Fe K$`\alpha `$ line, which we attribute to the continuum fit used by Brinkmann et al. (1998). More specifically we find that the line is several times broader and about 70% stronger than what Brinkmann et al. (1998) found. Comparing the parameters describing the continuum, we see that the Compton “reflection” component in the ASCA spectrum of Brinkmann et al. is poorly constrained and has an extremely high value that implies the reprocessor subtends an area larger than 3 times the area of the sky. Although this could be interpreted as an indication that the primary X-ray source does not emit isotropically, it is also possible that the strength of the Compton “reflection” component was overestimated.

Such poor constraints on the Compton “reflection” component are not surprising, since the ASCA sensitivity drops dramatically at $`E>8`$ keV, making it practically impossible to measure this component. Thus, the continuum around the Fe K$`\alpha `$ line is elevated, taking away the broad wings of the line and reducing its overall strength. In contrast, in our BeppoSAX spectrum we are able to determine the shape of the high-energy continuum unambiguously with the help of the PDS. We find a considerably weaker Compton “reflection” strength than Brinkmann et al. (1998) and a correspondingly broader and stronger Fe K$`\alpha `$ line.

In view of the strength of the Compton “reflection” component measured from the BeppoSAX spectrum, we re-analyzed the the ASCA spectra to check for consistency. We extracted and fitted the spectra with the same methodology described in §2 and §3, keeping the Compton “reflection” component fixed to a strength corresponding to $`\mathrm{\Omega }=2\pi `$, as measured from the BeppoSAX spectrum. We measured the EW of the Fe K$`\alpha `$ line for fixed rest energies of 6.4 and 6.7 keV. We found that the measurements from the ASCA spectra are consistent with those from the BeppoSAX spectrum. In particular, we measure a spectral index of $`1.87_{0.07}^{+0.08}`$, while the Fe K$`\alpha `$ rest EW is poorly constrained and lies in the range 130–700 eV, as listed in Table 3.

## 4 Summary and Discussion

We have presented new BeppoSAX and ASCA observations of three luminous ($`L_{210\mathrm{keV}}10^{45}\mathrm{erg}\mathrm{s}^1`$) radio-loud AGN, as part of our ongoing systematic study of the X-ray properties of such objects. We find that in all three cases the main underlying continuum can be described as a power law with photon index $`\mathrm{\Gamma }1.85`$, plus excess absorption. In two sources (PKS 2349–01 and 4C 74.26) the Fe K$`\alpha `$ emission line is detected with rest-frame $`\mathrm{EW}`$ 200 eV but with large uncertainties, of a factor of 2; in 4C 74.26 the line is unresolved while in PKS 2349–01 the line properties depend on the model adopted for the continuum. In 4C 74.26, a strong ($`\mathrm{\Omega }2\pi `$) Compton “reflection” component is also detected with the PDS instrument on BeppoSAX. In PKS 2349–01, on the other hand, the data allow the presence of a Compton “reflection” component but cannot provide meaningful constraints on its strength. This is unfortunate because the strength of this component affects the inferred properties of the Fe K$`\alpha `$ line dramatically, as we explained in detail in §3.2.

The average photon index for the three sources studied here is $`\mathrm{\Gamma }=1.86`$ with a small dispersion of $`\sigma __\mathrm{\Gamma }=0.01`$. The average slope of the X-ray continuum is consistent with the average photon index previously measured with ASCA for a small sample of 7 QSRs of comparable luminosity, which had $`\mathrm{\Gamma }=1.80`$ and $`\sigma __\mathrm{\Gamma }=0.15`$ (, Sambruna et al.1999). This is also consistent with the average slope for a sample of radio-quiet quasars from ASCA measurements, $`\mathrm{\Gamma }=1.86`$ and dispersion $`\sigma __\mathrm{\Gamma }=0.20`$ (see Table 5 of , Sambruna et al.1999). The similarity of slopes for the radio-loud and radio-quiet sources confirms that no contribution from the beamed emission of an unresolved jet is present in the radio-loud QSRs of our study (as expected based on their radio morphology, see above). It also suggests a similar continuum emission process in the two classes of AGN. This result will need to be confirmed using larger samples of both radio-loud and radio-quiet sources. Moreover, a wider spectral coverage will also be needed in order to constrain the strength of the Compton “reflection” hump and disentangle its effects from the slope of the intrinsic continuum.

Another result of this paper is the excess absorption found in all three objects. In 4C 74.26, where the observed spectra extend to low energies, we have searched for the signature of a warm absorber, but have not found conclusive evidence for it, nor were we able to rule it out. There are two possible origins for the excess absorption. The first is local excess absorption in the Galaxy. As 4C 74.26 lies at relatively low latitudes ($`b=19^{}.5`$), it is conceivable that the extra absorption is associated with a Galactic molecular cloud. Alternatively, the absorber could be intrinsic to the source. In PKS 2349–01 and 3C 323.1, the absorber is most likely intrinsic to the source because of their high Galactic latitudes. This is in agreement with earlier findings (, Sambruna et al.1999) that excess neutral absorption is common in radio-loud AGN and the absorbing column does not appear to be correlated with the X-ray luminosity of the source (these three objects fall in the general area occupied by radio-loud AGNs in the $`N_\mathrm{H}vsL_\mathrm{X}`$ plot shown in Figure 6 of , Sambruna et al.1999).

We have also detected an Fe K$`\alpha `$ emission line in two (possibly three) sources of our sample. In order to compare the strengths of these lines with what is found in other radio-loud and radio-quiet AGNs we plotted the rest EWs against the 2–10 keV luminosity in Figure 6. This is an update of the figure presented in (Sambruna et al.1999) and Eracleous et al. (2000). The data points represent radio-loud AGNs (, Sambruna et al.1999; Eracleous et al., 2000; Eracleous, 2002, and this paper). The shaded area represents the “X-ray Baldwin effect” for radio-quiet Seyferts and QSOs and its 90% dispersion (Nandra et al., 1997c). While at $`L_\mathrm{X}<5\times 10^{44}\mathrm{erg}\mathrm{s}^1`$ radio-loud AGNs have systematically lower Fe K$`\alpha `$ EWs than radio-quiet AGNs of comparable luminosity, the situation at higher luminosities is quite unclear because of the large uncertainties in the measured EWs. There are two possible ways of viewing Figure 6: either (a) radio-loud AGNs follow their own “X-ray Baldwin effect,” which is systematically offset from that of radio-quiet AGNs, or (b) radio-loud AGNs do not show a Baldwin effect at all. Unfortunately, the large uncertainties in the measured EWs, especially at high luminosities, do not allow us to distinguish between these two possibilities. Better measurements of the Fe K$`\alpha `$ EWs of luminous radio-loud AGNs ($`L_\mathrm{X}>10^{45}\mathrm{erg}\mathrm{s}^1`$) are sorely needed.

The X-ray Baldwin effect for radio-quiet AGN was interpreted by Nandra et al. (1997c) as the result of a progressive ionization of the inner accretion disk as the X-ray luminosity increases. We speculate that the lack of a similar trend for radio-loud AGN indicates that the ionization structure of the disk of radio-loud AGNs does not change significantly as the X-ray luminosity increases. This could happen, for example, in a scenario where the dense, optically thick matter in the disk is photoionized by external illumination with a hard spectrum of low ionization parameter. If the ionization parameter in the Fe K$`\alpha `$ emitting region is extremely low, changes in the ionizing luminosity even by an order of magnitude are accompanied by a proportional increase in the line luminosity. This would be the case if the central accretion disk is an ADAF (or similar structure) which has a low radiative efficiency and produces an X-ray luminosity that is disproportionately low compared to the accretion rate and a relatively hard spectrum. Ionization of the geometrically thin disk exterior to the ADAF by hard X-rays from the ADAF itself would thus produce the desired effect.

The above picture is supported by a number of other observational results, namely the systematically lower Fe K$`\alpha `$ EWs and FWHMs and weaker Compton “reflection” humps of radio-loud AGNs compared to their radio-quiet counterparts (see Woźniak et al., 1998; Eracleous et al., 2000; Grandi, 2002). Both of these differences can be understood in the context of a scenario in which the inner accretion disk is an ADAF rather than a geometrically thin and optically thick flow. The difference in Fe K$`\alpha `$ EWs and Compton “reflection” strengths is then a geometrical effect resulting from the small solid angle subtended by the outer disk to the primary X-ray source (Chen & Halpern, 1989; Zdziarski, Lubiński, & Smith, 1999). The difference in Fe K$`\alpha `$ FWHM follows simply from Kepler’s law.

To make further progress on this problem, more and better hard X-ray spectra of radio-loud AGN are needed, especially at high luminosity. The ideal instrument for this job is XMM-Newton. These spectra will serve to test the X-ray Baldwin effect for this class of objects. An independent test of the above hypothesis will be afforded by the profiles of the Fe K$`\alpha `$ lines, which would be expected to be significantly different in width and asymmetry from those of Seyfert galaxies.

We thank the anonymous referee for very helpful comments. CAH was funded by a fellowship from the Schreyer Honors College of the Pennsylvania State University and by NASA grants NAG5-9133 and NAG5-27017. This work was also partially supported by NASA grants NAG5-8369 and NAG5-10817. We have made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.