text:
# Large-System Phase-Space Dimensionality Loss in Stationary Heat Flows

## 1 Introduction

Aoki and Kusnezov stressed the usefulness of the “$`\varphi ^4`$” model in studying the size-dependence of classical heat conductivity simulations . They used this model to show that the entropy production associated with steady heat flow becomes extensive, ($`N`$), as the system size $`N`$ increases. When such heat-flow simulations are driven with two time-reversible Nosé-Hoover thermostats  imposing two temperatures, “hot” and “cold”, the steady-state phase-space distribution function becomes a multifractal object . The information dimension of this multifractal object lies below the corresponding equilibrium one by an extensive ($`N`$) “dimensionality loss” $`\mathrm{\Delta }D`$. The dimensionality loss, like the closely-related thermodynamic entropy production $`\dot{S}`$, is approximately quadratic in the deviation from equilibrium brought about by the temperature gradient $`T`$:

$$\mathrm{\Delta }DN(\mathrm{ln}T)^2.$$

This quadratic dependence is the usual prediction of linear transport theory. Its evaluation quantifies the rarity of nonequilibrium stationary states.

About ten years ago  we studied phase-space dimensionality loss in a family of two-dimensional shear flows. Those simulations were limited to relatively small systems. In that shear-flow work we were able to find systems for which the dimensionality loss barely, but significantly, exceeded the total dimensionality associated with the coordinates, momenta, and friction coefficients $`\{q,p,\zeta \}`$ of the thermostated boundaries.

At that time we stated that our data definitely showed that nonequilibrium stationary states can exhibit a reduced dimensionality, relative to equilibrium, even when projected into the subspace occupied by purely-Hamiltonian degrees of freedom. This interpretation was controversial . Many researchers had difficulty in accepting that nonequilibrium stationary states are typically associated with reduced dimensionality. In order to confirm that our interpretation was correct, we recently studied dimensionality reduction for stationary flows using the (Hamiltonian) $`\varphi ^4`$ model. Because heat transfer using this model requires only the simplest of boundaries—a single hot degree of freedom and a single cold one, in the simplest one-dimensional case—it turned out to be easy to get relatively large dimensionality losses .

A simple but convincing two-dimensional case involves the heat transfer from a single hot particle (with five phase-space coordinates $`\{x_H,y_H,p_{xH},p_{yH},\zeta _H\}`$) to a single cold one (with its own five coordinates) through a medium of $`N2=23`$ Hamiltonian particles (with 92 more phase-space coordinates, for a total of $`4N+2=102`$ phase-space coordinates). Typical $`(x,y)`$ particle trajectories in such a 25-particle system are shown in Fig. 1. In the

case of $`N=16`$ particles discussed in Ref. 2 the phase-space strange attractor turned out to have a dimensionality loss of 12.5, relative to the equilibrium distribution, exceeding the number of thermostated phase-space coordinates (5 hot and 5 cold) by 2.5.

The present work has two goals, (i) characterizing the overall loss of phase-space dimensionality $`\mathrm{\Delta }D`$, and (ii) estimating the corresponding loss $`\mathrm{\Delta }D_{}`$ in the Hamiltonian unthermostated part of phase space. This two-part study requires the numerical evaluation of the instantaneous Lyapunov exponents $`\{\lambda (t)\}`$, as well as their associated vectors $`\{\delta \}`$. The required techniques are sketched in Sec. II. In Sec. III we describe the $`\varphi ^4`$ model. In Sec. IV we undertake a systematic study of the number-dependence of the phase-space dimensionality loss $`\mathrm{\Delta }D`$. We consider a family of systems of the type illustrated in Fig. 1. The numerical work described in Sec. IV shows convincingly a $`\sqrt{1/N}`$ deviation from the large-system limit, with an extrapolated dimensionality loss of as much as four times (40) the number of dimensions (10) associated with the thermal driving mechanism. Next, in Sec. V, we study the projection of this dimensionality loss into the purely Hamiltonian subspace. We find that nearly all of the loss occurs in that subspace despite the purely Hamiltonian form of the equations of motion there. We discuss the projection technique and the importance of rapid rotation in phase space  to an understanding of the strange-attractor’s dimensionality reduction. The conclusions which follow from this work take up the final section.

## 2 Lyapunov Exponent Calculations

The usual “Lyapunov exponents”, $`\{\lambda _j\}`$, are time averages, over a sufficiently long time for convergence, of “local” (instantaneous) exponents $`\{\lambda _j(t)\}`$,

$$\lambda _j\lambda _j(t).$$

Here, and in what follows, we use angular brackets $`\mathrm{}`$ to indicate long-time averages. Each of the instantaneous exponents has associated with it an “offset” vector $`\delta _j`$ in phase space which describes the direction in which the growth or decay of phase-space separation is measured. If the offset vectors are infinitesimal in length, as in the present work, they can be replaced by a parallel set of unit vectors in “tangent space”.

The time evolution of the offset vector directions is governed by a continuous Gram-Schmidt orthonormalization which forces the vectors to remain mutually perpendicular and to evolve at fixed length. The two types of constraints, orthogonality and fixed length, can be imposed by a triangular array of Lagrange multipliers, $`\mathrm{\Lambda }_{ij}`$, where $`1jin`$ in an $`n`$-dimensional phase space . When the usual Cartesian coordinates are used the orthonormal offset vetors rotate rapidly .

The usual Lyapunov exponents $`\{\lambda _j\}`$ are the long-time-averaged values of the diagonal Lagrange multipliers $`\{\mathrm{\Lambda }_{jj}\}`$ required to enforce the orthonormalization constraints:

$$\lambda _j=\lambda _j(t)\mathrm{\Lambda }_{jj}.$$

At any instant of time the instantaneous Lyapunov exponents $`\lambda (t)`$ represent the orthogonal measurements of the $`n`$ growth (or decay) rates determined by the dynamical matrix $`D`$. $`D`$, an $`n\times n`$ matrix, is itself a phase function. The $`i`$th row of $`D`$ is made up of the $`n`$ derivatives of the $`i`$th equation of motion with respect to the $`n`$ phase variables.

The eigenvalues and eigenvectors of the dynamical matrix are relatively complicated . The eigenvalues are mainly rapidly varying complex-conjugate pairs, with the imaginary parts frequently vanishing, at phase-space singularities corresponding to parallel eigenvectors. Despite this complexity the orthonormal basis provided by the offset vectors makes it possible to measure smooth and well-defined growth rates. The time reversibility of the equations of motion guarantees related reversibility properties for $`D`$ and its eigenvectors and eigenvalues. Nevertheless, the past-based Lyapunov vectors governed by $`D`$ show a time-symmetry breaking intimately related to the second law of thermodynamics .

Consider a simple textbook (Sec. 5.4) illustration of these ideas, the motion of a driven thermostated particle in one dimension. The motion takes place in the three-dimensional $`(q,p,\zeta )`$ phase space with three Lyapunov exponents:

$$\begin{array}{ccccc}\lambda _1& =& \lambda _1(t)& =& \mathrm{\Lambda }_{11}\\ \lambda _2& =& \lambda _2(t)& =& \mathrm{\Lambda }_{22}\\ \lambda _3& =& \lambda _3(t)& =& \mathrm{\Lambda }_{33}.\end{array}$$

The equations of motion (with unit mass, force constant, and relaxation time) are

$$\dot{q}=p;\dot{p}=+1\zeta p;\dot{\zeta }=p^21,$$

for which the long-time solution is an attractor:

$$\{q(t),p(t),\zeta (t)\}\{t,+1,+1\}.$$

The matrix of equation-of-motion derivatives $`D`$ is

$$D=\left(\begin{array}{ccc}0& +1& 0\\ 0& \zeta & p\\ 0& 2p& 0\end{array}\right).$$

On the attractor the momentum $`p`$ and friction coefficient $`\zeta `$ can be replaced by their limiting values, (+1,+1). The three tangent-space $`\delta `$ vectors follow the ordinary differential equations:

$$\begin{array}{cccccccc}\dot{\delta }_1=& D\delta _1& & \mathrm{\Lambda }_{11}\delta _1& & & & \\ \dot{\delta }_2=& D\delta _2& & \mathrm{\Lambda }_{21}\delta _1& & \mathrm{\Lambda }_{22}\delta _2& & \\ \dot{\delta }_3=& D\delta _3& & \mathrm{\Lambda }_{31}\delta _1& & \mathrm{\Lambda }_{32}\delta _2& & \mathrm{\Lambda }_{33}\delta _3,\end{array}$$

where the six Lagrange multipliers follow easily from the time derivatives of the six orthonormality conditions:

$$\delta _1^2=\delta _2^2=\delta _3^21;$$

$$\delta _1\delta _2=\delta _2\delta _3=\delta _3\delta _10.$$

The steady-state time variation of the Lagrange multipliers, and the three $`\delta `$ vectors are shown in Fig. 2. The three Lyapunov exponents, like their three

instantaneous values $`\{\lambda _j(t)=\mathrm{\Lambda }_{jj}\}`$, sum to $`1`$:

$$\frac{\dot{q}}{q}+\frac{\dot{p}}{p}+\frac{\dot{\zeta }}{\zeta }\lambda _1(t)+\lambda _2(t)+\lambda _3(t)=1.$$

Note that the largest (time-averaged) Lyapunov exponent is 0 in this case, corresponding to a phase-space displacement in the direction of the motion $`\delta _1=(1,0,0)`$. The remaining five Lagrange multipliers, as well as the vectors $`\delta _2`$ and $`\delta _3`$, oscillate periodically, with a period of $`2\pi \sqrt{4/7}`$.

An analogous computation for an $`n`$-dimensional set of vectors involves computational work of order $`n^3`$. Nonequilibrium steady states generate chaotic attractors rather than the simple fixed point of the damped oscillator example. Nevertheless, the basic steps are the same: (i) propagating a phase-space reference trajectory; (ii) simultaneously propagating offset vectors in the neighborhood of the reference, with Lagrange multipliers and/or Gram-Schmidt orthonormalization imposing orthonormality; and (iii) averaging the diagonal multipliers to find the Lyapunov spectrum. Let us now consider the many-body model for which we carry out such investigations.

## 3 $`\varphi ^4`$ Heat Flow in Two Dimensions

In the “$`\varphi ^4`$” model we use here  we choose a square nearest-neighbor harmonic lattice with a quadratic Hooke’s Law pair potential for neighboring particles $`i`$ and $`j`$,

$$\varphi (r)=\frac{\kappa _2}{2}(rd)^2;r=|r_ir_j|>0.$$

In addition, each particle is tethered to its lattice site with a quartic potential, $`\frac{\kappa _4}{4}\delta r^4`$. The quartic tethers have two nice consequences. First, they provide “external forces”, and so prevent momentum conservation and ballistic energy transport. Second, they can provide chaos, so that the dynamics can have one or more positive Lyapunov exponents. The two consequences together can give Fourier heat conduction, even in one or two space dimensions. The heat conductivity for the $`\varphi ^4`$ model remains finite in the large-system limit, unlike many models, for which the conductivity vanishes or diverges in the large-system limit.

Here we simulate nonequilibrium heat-conducting stationary states by imposing thermostating forces on two of the $`N`$ particles. The thermostating forces are Nosé-Hoover  feedback forces $`\{\zeta p\}`$, linear in both the time-reversible friction coefficients $`\{\zeta \}`$ and the momenta $`\{p\}=\{(p_x,p_y)\}`$. The thermostated equations of motion for the hot and cold particles (one of each) are:

$$\begin{array}{ccccccccc}\{& m\dot{x}& =& p_x\hfill & ;& m\dot{y}& =& p_y\hfill & \}\\ \{& \dot{p}_x& =& F_x\zeta p_x\hfill & ;& \dot{p}_y& =& F_y\zeta p_y\hfill & \}\end{array}$$

$$\dot{\zeta }_{(H\mathrm{or}C)}=[(p^2/2mkT_{(H\mathrm{or}C)})1]/\tau ^2,$$

where $`\tau `$ is the characteristic response time of the thermostat forces, $`\{\zeta p\}`$. The full phase space describing this $`N`$-body two-dimensional system has $`4N+2`$ dimensions, with the extra two corresponding to the hot and cold friction coefficients $`\zeta _H`$ and $`\zeta _C`$.

For convenience in numerical work we choose the particle mass $`m`$, the spring constants $`\kappa _2`$ and $`\kappa _4`$, and the nearest-neighbor lattice spacing $`d`$ all equal to unity. The remaining parameters to set are the hot and cold temperatures, which we arbitrarily choose equal to the values from Ref. 2,

$$kT_H=0.009;kT_C=0.001,$$

and the thermostat relaxation times. We vary these times in the numerical work, but with the simplifying restriction that $`\tau _H`$ and $`\tau _C`$ have a common value, which we denote as $`\tau `$. We use the classic fourth-order Runge-Kutta integrator throughout, with a timestep $`dt=0.001\mathrm{or}0.002`$. To avoid numerical errors we have generated and compared results from two fully-independent simulation codes, one written in Vienna and the other written in Livermore.

## 4 Numerical Results—$`\mathrm{\Delta }D`$

The numerical evaluation of the dimensionality loss $`\mathrm{\Delta }D`$ is based on the connection between the Lyapunov spectrum and the dimensionality of the phase-space strange attractor. The Lyapunov exponents give the time-averaged relative growth and decay rates of the principal axes of a comoving infinitesimal phase-space hypersphere (or “extension in phase”). Kaplan and Yorke conjectured, evidently correctly , that a partial sum of these exponents (beginning with the largest one and proceeding toward the most negative one) changes from positive to negative when the (linearly-interpolated) number of terms in the sum is equal to the dimensionality of the phase-space strange attractor. This conjecture is “almost obvious”. It is evident that the (hyper)volume of any phase-space object with a positive sum of time-averaged Lyapunov exponents must diverge. Likewise a negative sum of time-averaged exponents indicates a vanishing hypervolume at long times. Any stationary process must generate an attractor which neither vanishes nor diverges.

The main difficulty in computing Kaplan-Yorke information dimensions is the unfavorable time-dependence associated with constraining the phase-space offset vectors to remain perpendicular to one another. With Gram-Schmidt orthonormalization in an $`n`$-dimensional phase space $`n`$ vectors, with $`n`$ components each, must all be propagated in time for sufficiently long that the time-averaged growth rates have converged. The computational work in orthogonalizing $`n(n1)/2`$ pairs of $`n`$-dimensional vectors varies as the cube of the number of particles, so that present computer speeds and processor numbers allow us to follow no more than a few hundred particles. In the present work we consider the simplest possible square systems of from $`4(2\times 2)`$ to $`144(12\times 12)`$ particles. The relaxation time $`\tau `$ is a free parameter. We choose it in the range $`1\mathrm{}8`$. Results become insensitive to $`\tau `$ once $`\tau `$ exceeds 6, which is comparable to the inverse Debye frequency of the lattice. Representative results for the largest Lyapunov exponent $`\lambda _1`$ and the dimensionality loss from the Kaplan-Yorke conjecture, $`\mathrm{\Delta }D`$, are given in Table 1.

Our own previous work, on color conductivity and shear flow , strongly suggests deviations in the Lyapunov spectrum of order the inverse system width, $`\sqrt{1/N}`$ in two dimensions. The present results are roughly consistent with this finding though an even slower variation with $`N`$ provides a comparable fit. We know of no previous systematic study of the variation of the Lyapunov spectrum with $`\tau `$. A cursory investigation shows that the dimensionality loss varies roughly as $`\tau ^1`$. Thus the dimensionality loss, with fixed boundary temperatures of 0.009 and 0.001, can be represented by

$$\mathrm{\Delta }D=42.7(5)\frac{89(2)}{\sqrt{N}}\frac{13.2(5)}{\tau }.$$

The standard deviation affecting the last digits of the fit parameters are given in brackets. We compare the second expression with our numerical data in Fig. 3. The rather good fit

indicates that the dimensionality loss persists in the large-system limit. We investigate the loss further in the following section.

## 5 Hamiltonian Projection: $`\mathrm{\Delta }D_{}`$

Ten years ago  we studied the rotation rate of the $`n`$ phase-space offset vectors $`\{\delta _j\}`$. We found that the rotation rate increases very rapidly with system size, soon becoming very large relative to the Lyapunov exponents themselves. This observation suggests that the instantaneous growth and decay rates in phase space—the instantaneous Lyapunov exponents—might become isotropic in the large-system limit. This suggests that the measured growth and decay rates are also closely associated with the subspaces spanned by the corresponding phase-space offset vectors $`\{\delta \}`$. Because the instantaneous Lyapunov exponents measure radial expansion and contraction, without any explicit rotational contribution, the instantaneous growth rates $`\{\lambda _j(t)=\mathrm{\Lambda }_{jj}\}`$ associated with every one of the phase-space directions contributing to a particular vector are identical. The contributions in a fixed phase-space direction are the summed-up contributions from the entire set of Lyapunov vectors.

At a particular phase-space point these contributions of the principal axes of the dynamical matrix to the Lyapunov exponents are all proportional to $`\mathrm{cos}^2(\theta )`$, where $`\theta `$ is the angular difference between the principal axis and the corresponding instantaneous $`\delta `$ vector. Because the Lyapunov exponents measure growth or decay in the direction of $`\delta `$, the same logarithmic growth rate, $`(d/dt)\mathrm{ln}\delta `$, applies to each component of the vector.

In the full phase space the information dimension of the strange attractor is given by the number of Lyapunov exponents whose sum is zero:

$$\underset{j=1}{\overset{k}{}}\lambda _j=0D=k;\mathrm{\Delta }D=4N+2k.$$

It has to be emphasized that $`k`$ (as well as $`k_{}`$ introduced below) is not an integer. A linear interpolation between two successive values of the Lyapunov sum is implied, with $`k`$ chosen such that the interpolated sum is precisely zero. In the Hamiltonian subspace the Lyapunov exponents contribute according to their projections into that space:

$$\underset{j=1}{\overset{k_{}}{}}\mathrm{cos}^2(\theta _j)\lambda _j(t)=0D_{}=\underset{j=1}{\overset{k_{}}{}}\mathrm{cos}^2(\theta _j);$$

$$\mathrm{\Delta }D_{}=\underset{j=k_{}}{\overset{4N+2}{}}\mathrm{cos}^2(\theta _j)=4N8\underset{j=1}{\overset{k_{}}{}}\mathrm{cos}^2(\theta _j).$$

The $`\mathrm{cos}^2(\theta )`$ form of the projection is required by the condition that changing the sign of the offset vectors leaves the growth rate unchanged. This form follows from the quadratic form describing an infinitesimal phase-space hyperellipsoid centered on a moving trajectory point. Note also that this weighting satisfies the normalization of the projections:

$$\underset{j=1}{\overset{4N+2}{}}\mathrm{cos}^2(\theta _j)4N8.$$

Partial sums give an effective number of exponents in the Hamiltonian projection of the full phase space. Thus the analog of Kaplan and Yorke’s conjecture for the Hamiltonian subspace is the effective number of exponents, $`\mathrm{cos}^2(\theta )`$, at which the projected sum, $`\lambda (t)\mathrm{cos}^2(\theta )`$, vanishes.

We explored the Kaplan-Yorke analog for the present problem, computing the delta vectors and their projections. The results are interesting. First, we noted that the summed-up local Lyapunov exponents are not strongly correlated with the directions in the subspace:

$$\lambda \mathrm{cos}^2(\theta )\lambda (t)\mathrm{cos}^2(\theta ).$$

This suggests that all the time-averaged projections of the various vectors into the Hamiltonian subspace are similar,

$$\mathrm{cos}^2(\theta )(4N8)/(4N+2).$$

See Figs. 4 and 5.

Fig. 4 demonstrates that the projection of the vectors becomes increasingly uniform as system size is increased, and is quite close to the average value, $`(4N8)/(4N+2)`$ expected for fully random projection directions. Fig. 5 compares the two estimates for the projected Lyapunov sums,

$$\lambda \mathrm{cos}^2(\theta )\mathrm{and}\lambda (t)\mathrm{cos}^2(\theta ).$$

The two sums vanish at nearly the same projected dimensionality, $`\mathrm{cos}^2(\theta )`$, indicating that the correlation of the exponents with direction is small.

We have used the analog of the Kaplan-Yorke formula to estimate the dimensionality reduction in the Hamiltonian subspace, $`\mathrm{\Delta }D_{}`$, and show these results in Fig. 6. Just as in the full phase space, the loss

of dimensionality varies smoothly with relaxation time and system size, and may be represented by

$$\mathrm{\Delta }D_{}=41.2(5)\frac{94(3)}{\sqrt{N}}\frac{10.1(5)}{\tau }.$$

It is worth pointing out that a naïve approach to dimensionality loss in the Hamiltonian subspace could be based on an orthonormalization of the Hamiltonian subspace only, propagating $`\delta `$ vectors with the Hamiltonian equations of motion but using the thermostated equations of motion for the underlying reference trajectory. This approach, which we explored years ago and which has recently been considered independently by Ken Aoki \[private communications, 2002\], is equivalent to considering the Lyapunov spectrum for a Hamiltonian system subject to time-dependent forces $`\{F(t)\}`$, where the forces are due to the thermostated particles. Such equations of motion, being Hamiltonian, satisfy Liouville’s Theorem , and guarantee that the corresponding Lyapunov spectrum is made up of pairs of positive and negative exponents, with zero sum.

## 6 Conclusions

The “$`\varphi ^4`$” model shows conclusively that the nonequilibrium steady-state loss of phase-space dimensionality can easily exceed the dimensionality associated with the system boundaries. The present results also confirm that the deviations from the large-system limit vary according to a power law, $`N^p`$, $`1/4p1/2`$, roughly compatible with the inverse of the system size, $`1/\sqrt{N}`$ in two dimensions. For $`\tau 5`$ they have in addition deviations in the relaxation frequency $`(1/\tau )`$. The quartic tethers are a particularly useful feature of the model, which make it possible to carry out the simulations without the need to take thermal expansion explicitly into account.

The dimensionality loss $`\mathrm{\Delta }D`$ is expected to be extensive as has been stressed in the introduction. For fixed thermostat temperatures $`T_H`$ and $`T_C`$ the temperature gradient is determined by the system size, proportional to $`1/\sqrt{N}`$ in two dimensions, and $`\mathrm{\Delta }D`$ becomes independent of $`N`$. This is indeed the case in the large-particle limit $`N\mathrm{}`$. Our simulations thus confirm the extensivity of the dimensionality reduction for stationary heat flow in the linear-response limit. Far from equilibrium weak deviations from this $`N`$-dependence are found.

Table 1 reveals that for short Nosé-Hoover relaxation times the dimensionality loss $`\mathrm{\Delta }D_{}`$ in the purely Hamiltonian part of phase space may even slightly exceed the dimensionality loss $`\mathrm{\Delta }D`$ in the full phase space. This is because $`\mathrm{\Delta }D_{}`$ is overestimated by one due to the lack of energy conservation in the Hamiltonian subspace once the thermostats are added.

The projection of the phase-space offset vectors $`\{\delta \}`$ into the Hamiltonian subspace developed here shows that most of the dimensionality loss occurs in a part of the system which obeys purely conservative Hamiltonian equations of motion. Rapid rotation is responsible. This rotation nearly eliminates the correlations between phase-space contraction and direction. Evidently phase-space contraction is not only real, but relatively simple, and certainly must persist in the large-system thermodynamic limit. Thus the present results corroborate our interpretation of the second law of thermodynamics for nonequilibrium stationary states . Such states occupy not just a reduced volume in phase space. They are restricted to a subspace of reduced dimensionality, with the dimensionality loss simply related to the rate of external thermodynamic entropy production. Useful models illustrating dimensionality loss for periodic color conducting or shear flows can now be developed as extensions of this idea.

## Acknowledgements

Work at the University of Vienna was supported by the Austrian Fonds zur Förderung der wissenschaftlichen Forschung, Project P-15348. Work at Livermore was performed under the auspices of the United States Department of Energy through University of California Contract W-7405-Eng-48. We very much appreciate the interest and encouragement of Kenichiro Aoki, Aurel Bulgac, Carol Hoover, Dimitri Kusnezov, and David Ruelle.

w


format: text

text:
# Coexistence of Coherence and Incoherence in Nonlocally Coupled Phase Oscillators: A Soluble Case

## 1 Introduction

Large populations and continuous fields of coupled oscillators form a representative class of synergetic systems met in a wide range of scientific disciplines from physics, chemistry, engineering, biology to brain science . Collective dynamics of coupled oscillators depends crucially on the range of their mutual coupling. It has recently been realized that when the coupling is nonlocal, the patterns which emerge could be drastically different from those which we expect for oscillators with local or global coupling . The implication of this fact is relevant even to what we conventionally call locally coupled systems, typically reaction-diffusion systems. This is because it may happen that nonlocality can arise effectively as a result of elimination of some variables, e.g., rapidly diffusing components in the case of reaction-diffusion dynamics. Among the variety of patterns which are characteristic to nonlocally coupled oscillators, we will focus our attention below on a particular class of patterns in which the whole medium is separated into two domains of qualitatively different dynamics. Specifically, the oscillators are mutually synchronized in one domain while they are completely desynchronized. in the other domain. A preliminary work on such dynamics was reported recently. We will present below a more thorough investigation of this problem.

The collective dynamics of our concern is similar to the collective synchronization in globally coupled oscillators with distributed natural frequencies  where the whole population splits into two subpopulations each composed of synchronized and desynchronized oscillators. There the systems is stationary in a statistical sense within a constant drift of the collective phase corresponding to the oscillation of the population as a whole. Similarly to the case of global coupling, the individual oscillators can be regarded as being controlled by a common mean field, although the latter is now space dependent. The oscillators under consideration are identical in nature, still we may develop an exact mean-field theory similar to the one developed for globally coupled oscillators. Commonly to these theories, a suitably defined order parameter representing the collective state is determined from a condition for self-consistency which must exist between the collective dynamics and the dynamics of the individual oscillators.

In section 2, some results of our numerical simulation will be presented. The simulation is carried out first on a nonlocal version of the complex Ginzburg-Landau (CGL) equation in one space dimension with periodic boundary conditions, and then on another equation obtained through the phase reduction of the nonlocal CGL. In section 3, our theory will be developed and a functional self-consistency equation will be derived for the order parameter. Then the numerical solution of the self-consistency equation will be compared with the simulation results presented in Section 2. Concluding remarks will be given in the final section.

## 2 Coexistence of coherence and incoherence

As a simplest model for densely and uniformly distributed oscillators with nonlocal coupling, let us consider the following equation for a complex amplitude $`A`$ which we call nonlocally coupled complex Ginzburg-Landau equation. In one spatial dimension our model is expressed as

$$\frac{}{t}A(x,t)=(1+i\omega _0)A(1+ib)|A|^2A+K(1+ia)\left(Z(x,t)A(x,t)\right).$$
(1)
The quantity $`Z(x,t)`$, which we call the mean field, represents the effects of the nonlocal coupling and is given by

$$Z(x,t)=G(xx^{})A(x^{},t)𝑑x^{}.$$
(2)
The coupling function $`G`$ changes with the distance as

$$G(y)=\frac{\kappa }{2}\mathrm{exp}(\kappa |y|)$$
(3)
and is normalized. The exponential form of $`G`$ is in fact a natural consequence of the reductive derivation of Eq. (1) from a certain class of reaction-diffusion systems when the latter involve an inactive diffusive component to be eliminated adiabatically. Note that the last term in Eq. (1) representing the coupling is so arranged that it may vanish when $`A`$ is uniform in space. If $`A`$ is nonuniform but sufficiently long-waved, this term can be approximated with a diffusion term, and then Eq. (1) is reduced to the ordinary CGL.

Equation (1) is numerically simulated on a finite interval $`[0,1]`$ with periodic boundary conditions. In a suitable range of parameters, in which the uniform oscillation is linearly stable, the system can develop peculiar patterns when the initial perturbation is finite. Figure 1 shows a typical example of such phase patterns. The initial condition is such that the modulus of $`A`$ equals 1 everywhere and its phase changing randomly still its envelope being nearly symmetric about the midpoint $`x=1/2`$. Since the coupling constant $`K`$ is assumed to be relatively small, the local oscillations attain almost full amplitude, so that we will mainly be interested in the phase pattern. The pattern consists of two distinct domains. In one domain which appears near the boundaries, the pattern is spatially continuous and smooth, whereas in the central domain spatial continuity seems completely lost. The pattern as a whole looks steadily advancing upward, implying that the entire population is undergoing regular collective oscillation with a definite frequency $`\mathrm{\Omega }`$.

To investigate the nature and origin of such patterns in further detail, Eq. (1) is reduced to a phase equation which would be much easier to analyze. The phase reduction is valid when the coupling strength $`K`$ is small, which is actually the case, and the reduced equation takes the form

$$\frac{}{t}\varphi (x,t)=\omega G(xx^{})\mathrm{sin}\left(\varphi (x,t)\varphi (x^{},t)+\alpha \right)𝑑x^{},$$
(4)
where the time scale has been so changed as to normalize the coupling strength, $`\omega `$ is the natural frequency rescaled accordingly, and the phase constant $`\alpha `$ in the phase-coupling function is related to the original parameters through

$$\mathrm{tan}\alpha =\frac{ba}{1+ab},\alpha (ba)>0.$$
(5)
Under the parameter condition corresponding to the previous numerical simulation of Eq. (1), a phase pattern similar to FIG. 1 was obtained, which is displayed in FIG. 2.

We now focus on the phase patterns obtained for our phase equation. Let us introduce relative phase $`\psi `$ by

$$\varphi =\mathrm{\Omega }t+\psi $$
(6)
which describes the dynamics of the phase deviation from the reference motion with some drift velocity $`\mathrm{\Omega }`$ whose value is still open. We rewrite Eq. (4) using $`\psi `$ as

$$\frac{}{t}\psi (x,t)=\omega \mathrm{\Omega }𝑑x^{}G(xx^{})\mathrm{sin}\left(\psi (x,t)\psi (x^{},t)+\alpha \right).$$
(7)
As a generalization of the theory of synchronization transition of globally coupled phase oscillators with frequency distribution, we introduce a complex order parameter with modulus $`R`$ and phase $`\mathrm{\Theta }`$ through

$$𝑑x^{}G(xx^{})\mathrm{exp}[i\psi (x^{},t)]=R(x,t)\mathrm{exp}[i\mathrm{\Theta }(x,t)].$$
(8)
Unlike the order parameter defined for globally coupled oscillators, the above quantity is space-dependent. Still a physical picture is valid such that we are practically working with an assembly of independent oscillators under the control of a common forcing field represented by $`R`$ and $`\mathrm{\Theta }`$. This is confirmed by rewriting Eq. (7) in terms of the order parameter into a forced one-oscillator equation as

$$\frac{}{t}\psi (x,t)=\omega \mathrm{\Omega }R(x,t)\mathrm{sin}\left(\psi (x,t)+\alpha +\mathrm{\Theta }(x,t)\right).$$
(9)
The spatial profiles of $`R(x)`$ and $`\mathrm{\Theta }(x)`$ after a long-time average are displayed in FIG. 3a and 3b, respectively. We see that the forcing mean-field amplitude is stronger near the boundaries and weaker near the center of the system. If the phase distribution such as shown in FIG. 2 is statistically stationary within a constant overall drift with velocity $`\mathrm{\Omega }`$, then $`R`$ and $`\mathrm{\Theta }`$ are time-independent. Their approximate constancy is confirmed by the fact that the fluctuation of the time sequence of $`R(x,t)`$ observed at the midpoint of the system $`x=1/2`$ gives a standard deviation as small as 0.008 which comes possibly from a finite-size effects. Similar property can also be confirmed for $`\mathrm{\Theta }(x,t)`$.

Figure 3c shows a distribution of the actual frequencies $`\overline{\omega }(x)`$ of the individual oscillators. This was also obtained through a long-time average. In the coherent domain, the oscillation frequencies have an identical value $`\mathrm{\Omega }`$, while in the incoherent domain they are distributed but give a well-defined continuous curve.

The dynamics which is going on is now clear. The system is divided into two subgroups of oscillators. In the first group, the forcing amplitude is large enough for the oscillators to be entrained, so that they oscillate with an identical frequency $`\mathrm{\Omega }`$. In the second group, in contrast, the forcing amplitude is too weak for entrainment, so that the frequencies of the individual oscillators differ from $`\mathrm{\Omega }`$. Since the actual frequencies of these desynchronized oscillators should depend on the local amplitude of the forcing, the nonuniform but smooth pattern of $`R(x)`$ as shown in FIG. 3a implies a smooth distribution of the actual frequencies, and this is consistent with FIG. 3c showing $`\overline{\omega }(x)`$. Our next problem is how to predict stationary patterns of $`R(x)`$ and $`\mathrm{\Theta }(x)`$ theoretically, and this will be the subject of the next section.

## 3 Theory

The basic equations to work with below are Eqs. (8) and (9). Assuming $`R`$ and $`\mathrm{\Theta }`$ to be time-independent, we first try to find solutions of Eq. (9) for the phases $`\psi (x,t)`$ each of which should be a function of $`R(x)`$ and $`\mathrm{\Theta }(x)`$. These solutions are then substituted into Eq. (8), leading to an equation for determining $`R`$ and $`\mathrm{\Theta }`$ as a functional of $`R`$ and $`\mathrm{\Theta }`$ in a self-consistent manner. Since the collective frequency $`\mathrm{\Omega }`$ is still to be determined, we will be working with a nonlinear eigenvalue problem.

Here a question may arise as to the seeming contradiction between the assumed time-independence of $`R`$ and $`\mathrm{\Theta }`$ and the apparent time-dependence of the solution of Eq. (9) which is true if $`R(x)`$ is below a certain critical value. The same problem was encountered in the problem of synchronization transition in globally coupled phase oscillators with frequency distribution. The answer to this question is that the factor $`\mathrm{exp}[i\psi (x^{},t)]`$ appearing in the integral of Eq. (8) should be replaced with is statistical average $`\mathrm{exp}[i\psi (x^{},t)]`$. Such an average can be calculated using a suitable invariant measure which is easy to find.

The successive steps to be taken for reaching the final solution are the following. Note first that Eq. (9) admits two types of solutions depending on the location of the oscillators of concern. They are constant solutions and drifting solutions. If $`|R/(\omega \mathrm{\Omega })|1`$, we have a time-independent solution given by

$$\theta _0(R)=\text{Sin}^1\left(\frac{\omega \mathrm{\Omega }}{R(x)}\right)\alpha .$$
(10)
If $`|R/(\omega \mathrm{\Omega })|>1`$, in contrast, we have a drifting solution. The actual frequency $`\overline{\omega }`$ of the corresponding desynchronized oscillator as a function or $`R`$ is given by

$$\overline{\omega }=\mathrm{\Omega }+(\omega \mathrm{\Omega })\sqrt{1\left(\frac{R}{\omega \mathrm{\Omega }}\right)^2},$$
(11)
which equals the collective frequency $`\mathrm{\Omega }`$ plus the mean drift velocity of $`\psi `$. Instead of using such drifting solutions in Eq. (9), we use an invariant measure, i.e., the probability density $`p(\theta ,R)`$ such that $`\psi `$ takes on value $`\theta `$ for given $`R`$, which is inversely proportional to the drift velocity of $`\psi `$ at $`\psi =\theta `$. Thus,

$$p(\theta ,R)=C(\omega \mathrm{\Omega }R\mathrm{sin}\theta )^1,$$
(12)
where C is a normalization constant given by $`C=(2\pi )^1\sqrt{(\omega \mathrm{\Omega })^2R^2}`$. As noted earlier, this invariant measure is used for replacing the exponential factor appearing on the left-hand side of Eq. (8) with its statistical average. The physical idea underlying this replacement is that an equilibrium statistical ensemble may be applicable to each small local subsystem because infinitely many oscillators which are independently oscillating are already contained in such a subsystem.

The spatial profile of $`R(x)`$ in FIG. 3 implies that there are spatial points $`x=\frac{1}{2}\pm x_c`$ at which $`R`$ takes a critical value separating the two types of solutions described above. Namely, in the inner domain ($`|x^{}\frac{1}{2}|x_c`$) the solutions are drifting, while in the outer domain ($`|x^{}\frac{1}{2}|x_c`$) they are constant in time. Thus Eq. (8) becomes

$$_0^1𝑑x^{}G(|xx^{}|)\mathrm{exp}\left[i\left(\mathrm{\Theta }(x^{})\alpha \right)\right]h\left(R(x^{})\right)=R(x)\mathrm{exp}[i\mathrm{\Theta }(x)],$$
(13)
where

$$h(R)=\{\begin{array}{cc}\mathrm{exp}[i\theta _0(R)]\hfill & (|x^{}\frac{1}{2}|x_c)\hfill \\ _0^{2\pi }\mathrm{exp}(i\theta )p(\theta ,R)𝑑\theta \hfill & (|x^{}\frac{1}{2}|<x_c)\hfill \end{array}$$
(14)
The above expression for $`h(R)`$ in each domain can further be made explicit as

$`\mathrm{exp}[i\theta _0(R)]`$ $`=`$ $`\sqrt{1\left({\displaystyle \frac{\omega \mathrm{\Omega }}{R}}\right)^2}+i{\displaystyle \frac{\omega \mathrm{\Omega }}{R}},`$ (15)
$`{\displaystyle _0^{2\pi }}\mathrm{exp}(i\theta )p(\theta ,R)\theta `$ $`=`$ $`{\displaystyle \frac{i}{R}}\left\{\omega \mathrm{\Omega }\sqrt{(\omega \mathrm{\Omega })^2R^2}\right\}.`$ (16)
Equations (13)$``$(16) constitute a functional self-consistency condition from which $`R(x)`$, $`\mathrm{\Theta }(x)`$ and $`\mathrm{\Omega }`$ are to be determined. By using $`R(x)`$ and $`\mathrm{\Omega }`$ thus obtained, the actual frequencies of the desynchronized oscillators are determined from Eq. (11).

There is a unique value $`\mathrm{\Omega }`$ of the collective frequency for which our functional self-consistency equation admits a solution. Such solution can be obtained numerically, and is compared with our numerical simulation. Some of the results are given in FIG. 4a$``$4c. We find that the agreement between the theory and numerical experiments is almost perfect for each quantity of $`R(x)`$, $`\mathrm{\Theta }(x)`$, $`\mathrm{\Omega }`$ and $`\overline{\omega }(x)`$.

## 4 concluding remarks

The present theory makes full use of the fact that the number of oscillators contained in the range of coupling is practically infinite. This fact allows us to apply a mean field theory, and our problem is reduced to a one-oscillator problem supplemented with a self-consistency condition. It is well known that the same idea works nicely for globally coupled oscillators. The present theory suggests that, beyond a special problem discussed in the present paper, the mean-field idea may work more generally for nonlocally coupled systems as far as the coupling radius remains much longer than the minimum length scale associated with the discreteness of the spatial distribution of the constituent oscillators.

Acknowledgement

The present work gives a slight extension of the my old work on globally coupled phase oscillators which was developed in the middle of 1970’s, shortly after which I spent a very fruitful year at Professor Hermann Haken’s Institute in Stuttgart. Therefore, it is my great pleasure to dedicate this small work to Professor Haken.

Figure Captions

* Instantaneous spatial distribution of the phases obtained from Eq. (1) approximated with an array of 512 oscillators distributed over an interval of length 1 with periodic boundary conditions. Parameter values are $`a=1.0`$, $`b=0.88`$, $`\kappa =4.0`$.
* Instantaneous spatial distribution of the phases obtained from Eq. (4) approximated with an array of 512 oscillators distributed over an interval of length 1 with periodic boundary conditions. Parameter values are $`\alpha =1.457`$, $`\kappa =4.0`$.
* Spatial profiles of the long-time average of the order parameter amplitude $`R(x)`$ (a), the order parameter phase $`\mathrm{\Theta }(x)`$ (b) and actual frequencies $`\mathrm{\Omega }(x)`$ (c), each obtained from Eq. (4).
* Theoretical curves for $`R(x)`$ (solid line in (a)), $`\mathrm{\Theta }(x)`$ (solid line in (b)) and $`\mathrm{\Omega }(x)`$ (broken line in (c)) are compared with the corresponding quantities in FIG. 3 obtained numerically.


text:
# Untitled Document

ANÁLISE TERMODINÂMICA DE UM CIRCUITO RLC

Rodrigo de Abreu

Centro de Electrodinâmica e Departamento de Física do IST

Abstract
We analyse a RLC circuit taking the Second Law of Thermodynamics into
consideration.

Introdução

Hoje em dia, torna-se cada vez mais difícil definir onde se situa a fronteira entre a Mecânica e a Termodinâmica. Em primeiro lugar, torna-se difícil em inúmeros problemas fazer uma distinção clara, e com significado físico, entre os conceitos de trabalho e calor associados à Primeira Lei da Termodinâmica . Em segundo lugar, como é sabido, uma das vias para a interpretação física de alguns conceitos termodinâmicos, como o calor e a entropia, é a partir da Mecânica Estatística. Por último, na maior parte dos problemas de Mecânica envolvendo sistemas dissipativos, a análise realizada debruça-se exclusivamente sobre os aspectos mecânicos, não sendo feita qualquer tentativa para uma abordagem termodinâmica do problema.

O problema que iremos aqui abordar é, na sua essência, extremamente simples, sendo descrito pela Lei de Hooke. Esta lei postula a proporcionalidade entre a força aplicada a uma mola e o elongamento sofrido por esta

$$F=k(L_1L_0),$$
(1)
onde $`k`$ é a constante de restituição da mola e $`L_0`$ e $`L_1`$ representam, respectivamente, os comprimentos inicial e final da mola. Uma mola que obedeça a (1), com $`k`$ constante em qualquer circunstância, é chamada aqui mola ideal. Uma outra situação mais real, onde a proporcionalidade anterior já não se verifique, irá ser tratada igualmente neste trabalho numa secção à parte.

1 $``$ A Mola Ideal como Sistema Dissipativo

Considere-se então em primeiro lugar uma mola ideal, de constante $`k`$ e comprimento $`L_0`$, à qual se aplica uma força constante $`F_0`$, por exemplo, pendurando-se na vertical uma massa de peso $`F_0`$=$`mg`$. Como é sabido, a mola sofre uma elongação até à sua nova posição de equilíbrio $`L_1`$, mesmo que a mola se encontre no vácuo, sendo este ponto simplesmente determinado por

$$k(L_1L_0)=F_0.$$
(2)
O trabalho realizado pela força aplicada exterior na elongação da mola é usando (2) dado por

$$\tau =_{L_0}^{L_1}F_0𝑑x=F_0(L_1L_0)=\frac{F_0^{\mathrm{\hspace{0.33em}2}}}{k},$$
(3)
enquanto que a energia potencial por esta adquirida é dada por

$$V=_{L_0}^{L_1}kx𝑑x=\frac{1}{2}k(L_1L_0)^2=\frac{1}{2}\frac{F_0^{\mathrm{\hspace{0.33em}2}}}{k},$$
(4)
isto é, igual a metade do trabalho $`\tau `$. Desta diferença emerge imediatamente a seguinte pergunta: Para onde vai então a parte restante do trabalho realizado pela força $`F_0`$? Não fica concerteza armazenado no sistema sob a forma de energia cinética, dado que a mola fica parada no fim do movimento.

A resposta a esta questão, reside obviamente no facto de termos considerado a elongação da mola como um processo sem atrito. Na realidade, a elongação da mola ocorre sempre com dissipação de energia e quaisquer que sejam as características da mola ou do meio circundante que a envolva, a energia dissipada é exactamente igual a metade do trabalho realizado pela força aplicada $`F_0`$. Isto é, mesmo que o meio circundante seja o vácuo, onde não existe atrito com o ar, a energia dissipada apresenta sempre o mesmo valor, $`W_{diss.}`$=$`\tau `$/2.

Uma situação de certo modo análoga surge quando se carrega um condensador de capacidade $`C`$, inicialmente descarregado, ligando-o a uma bateria de força electromotriz $`ϵ`$. Vamos supor que a carga do condensador é feita ligando-se em série, com o condensador e a bateria, uma resistência $`R`$, a qual inclui a resistência da própria bateria. A solução deste problema é bem conhecida dos cursos de Electromagnetismo. A diferença de potencial (d.d.p.) aos terminais do condensador é dada por

$$v_c(t)=ϵ(1e^{t/RC}),$$
(5)
pelo que o condensador vai acabar por ficar carregado com uma d.d.p. $`v_c(t)`$=$`ϵ`$, quando $`t`$$``$$`\mathrm{}`$, e portanto com uma energia electrostática final igual a

$$W_e=\frac{1}{2}Cϵ^2.$$
(6)
Por outro lado, a corrente que percorre o circuito, enquanto ele estiver a ser carregado, é dada por

$$i=\frac{ϵ}{R}e^{t/RC},$$
(7)
o que permite calcular imediatamente, tanto a energia dissipada por efeito de Joule na resistência $`R`$

$$W_{diss.}=_0^{\mathrm{}}Ri^2𝑑t=\frac{1}{2}Cϵ^2,$$
(8)
como o trabalho fornecido pela bateria

$$\tau =_0^{\mathrm{}}ϵi𝑑t=Cϵ^2.$$
(9)
Tem-se aqui também, tal como no caso da mola, que a energia dissipada na resistência é exactamente igual a metade do trabalho fornecido pela bateria; e este resultado é tanto mais notável quando se observa que o resultado a que se chega é independente do valor da resistência $`R`$. Isto é, se fizermos tender $`R`$$``$0, o condensador é carregado num tempo extremamente curto, com um valor inicial de corrente $`i(0)`$$``$$`\mathrm{}`$. Contudo, a energia total dissipada mantem-se sempre constante, apresentando o valor $`W_{diss.}`$=$`\tau `$/2.

Ora, a situação que aqui pretendemos tratar, de uma mola ideal sujeita a uma força constante $`F_0`$, apresenta semelhanças óbvias com o caso do condensador. Em primeiro lugar, em vez de nos preocuparmos em determinar apenas o ponto de equilíbrio dado por (2), temos de considerar que a equação que descreve a elongação da mola é na realidade dada por

$$m\frac{d^2x}{dt^2}+\beta \frac{dx}{dt}+kx=F_0,$$
(10)
sendo $`x`$=($`L`$$``$$`L_0`$) o valor da elongação da mola em cada instante e onde fenomenologicamente se introduziu uma força de atrito $`F_{at.}`$=$`\beta `$$`v`$, tendo em conta os efeitos dissipativos que ocorrem na mola. A equação (10) apresenta como solução, no caso de um regime oscilante amortecido

$$x(t)=\frac{F_0}{k}+Ae^{\lambda t}\mathrm{cos}(\omega t+\varphi ),$$
(11)
com $`\lambda `$=$`\beta `$/$`2m`$ e $`\omega `$=$`\sqrt{(k/m)\lambda ^2}`$. No caso da mola se encontrar parada no início, as condições de continuidade em $`t`$=0 para as energias potencial e cinética, determinam as constantes de integração: $`A`$=$``$($`F_0`$/$`k`$)$`\mathrm{sec}\varphi `$; $`\varphi `$=$``$$`\mathrm{arctan}(\lambda /\omega )`$. No limite $`t`$$``$$`\mathrm{}`$, tem-se como anteriormente ($`L_1`$$``$$`L_0`$)=$`F_0`$/$`k`$.

Os sistemas mecânico e eléctrico aqui descritos apresentam efectivamente o mesmo tipo de soluções. A aparente diferença entre as soluções (5) e (11) resulta simplesmente do facto de termos considerado nulo o coeficiente de indução do circuito. Se ao invés, tivessemos admitido a existência de um elemento indutivo no circuito, as soluções seriam também do tipo oscilante amortecido. Contudo, o balanço energético manter-se-ia inalterado, na medida em que o elemento indutivo não dissipa nem armazena energia no final. Note-se que após o condensador estar carregado $`i`$=0, pelo que a energia magnética no elemento indutivo é também nula.

Se multiplicarmos agora ambos os membros de (10) pela variação elementar $`dx`$ e integrarmos essa equação entre $`L_0`$ e $`L_1`$, obtemos a equação que traduz o balanço de energia na mola. O primeiro termo do membro de lado esquerdo é obviamente nulo, na medida que não existe variação de energia cinética entre os pontos de equilíbrio inicial e final, o segundo termo representa a energia dissipada na mola $`W_{diss.}`$, o terceiro termo representa a variação de energia potencial dada por (4), enquanto que o termo do lado direito é o trabalho (3) realizado pela força exterior. Este simples balanço permite-nos escrever para a energia dissipada no termo de atrito $`\beta `$

$$W_{diss.}=_{L0}^{L1}\beta \frac{dx}{dt}𝑑x=\frac{1}{2}F_0(L_1L_0).$$
(12)
É interessante notar ainda que a conclusão a que se chegou é independente do valor de $`\beta `$, podendo inclusivamente ter-se $`\beta `$$``$0. Este caso corresponde no circuito $`RC`$ à situação $`R`$$``$0, em que $`i(0)`$$``$$`\mathrm{}`$, mas onde a energia dissipada por efeito de Joule continua, mesmo neste limite, a ser dada por (8). Regressando agora ao caso da mola e supondo $`m`$$``$0, de forma ao sistema ser descrito igualmente por uma equação diferencial de primeira ordem (correspondente ao limite $`L`$$``$0 no circuito $`RLC`$), ao limite $`\beta `$$``$0 corresponde uma velocidade inicial infinita ($`dx`$/$`dt`$)0$``$$`\mathrm{}`$. Contudo, a energia dissipada na mola continua a ser dada por (12). Assim, embora este facto não seja normalmente tratado nos cursos elementares, uma mola ideal, sujeita a uma força aplicada $`F_0`$ constante, é sempre dissipativa. Este facto verifica-se independentemente das características materiais da mola.

2 $``$ Análise Termodinâmica da Mola Ideal

O resultado a que se chegou em (12) mostra que a energia dissipada na mola é sempre igual a metade do trabalho realizado pela força aplicada exterior, independentemente do valor do termo de atrito $`\beta `$. Assim, mesmo que o meio onde a mola se encontre seja o vácuo e que não exista atrito com o ar, há sempre um atrito interno da própria mola. A energia dissipada internamente dá origem a uma elevação de temperatura da mola acompanhada por uma transferência de energia para o ambiente, esta última sob a forma de radiamento no caso do meio envolvente ser o vácuo.

Este problema pode ser tratado agora sob o ponto de vista termodinâmico. Em primeiro lugar, a entropia da mola obedece à equação de estado $`S`$=$`S`$($`T`$,$`L`$), pelo que uma variação elementar desta é dada por

$$dS=\left(\frac{S}{T}\right)_LdT+\left(\frac{S}{L}\right)_TdL,$$
(13)
enquanto que para a energia interna da mola podemos escrever a relação termodinâmica fundamental para um processo infinitesimal

$$dU=TdS+FdL.$$
(14)
Como sabemos, a capacidade calorífica para um comprimento $`L`$ é dada por

$$C_L=\left(\frac{U}{T}\right)_L=T\left(\frac{S}{T}\right)_L,$$
(15)
o que nos permite obter para o primeiro termo de (13)

$$\left(\frac{S}{T}\right)_L=\frac{C_L}{T}.$$
(16)
Por outro lado, o segundo termo de (13) pode ser calculado usando uma das relações termodinâmicas de Maxwell, com as variáveis habituais pressão e volume substituídas aqui por $`F`$ e $`L`$, pelo que usando (1), onde $`k`$ é independente de $`T`$, podemos escrever

$$\left(\frac{S}{L}\right)_T=\left(\frac{F}{T}\right)_L=0.$$
(17)
Temos assim que neste caso a entropia pode ser calculada a partir da relação

$$dS=\frac{C_L}{T}dT,$$
(18)
a qual uma vez integrada, admitindo que $`C_L`$ não depende de $`T`$, permite obter para a variação de entropia, entre os estados inicial $`S`$=$`S`$($`T_0`$,$`L_0`$) e final $`S`$=$`S`$($`T_1`$,$`L_1`$), a seguinte expressão

$$\mathrm{\Delta }S=C_L\mathrm{ln}\left(\frac{T_1}{T_0}\right).$$
(19)
Neste caso, a transformação é isentrópica se $`T_1`$=$`T_0`$.

No que se refere agora à energia interna, usando (1) e (18) podemos escrever (14) sob a forma

$$dU=C_LdT+k(LL_0)dL,$$
(20)
obtendo para a variação de energia interna entre os mesmos estados inicial e final,

$$\mathrm{\Delta }U=C_L(T_1T_0)+\frac{1}{2}k(L_1L_0)^2.$$
(21)
Se admitirmos que não há transferência de energia para o ambiente, a variação de energia interna da mola $`\mathrm{\Delta }U`$ é igual ao trabalho realizado pela força aplicada exterior na elongação da mola (3)

$$\tau =\frac{F_0^2}{k},$$
(22)
sendo o comprimento final determinado por (2), ($`L_1`$$``$$`L_0`$)=$`F_0`$/$`k`$. O aumento de temperatura é dado assim por

$$T_1=T_0+\frac{1}{2}\frac{F_0^{\mathrm{\hspace{0.33em}2}}}{kC_L},$$
(23)
enquanto que para a variação de entropia tem-se

$$\mathrm{\Delta }S=C_L\mathrm{ln}\left(1+\frac{1}{2}\frac{F_0^{\mathrm{\hspace{0.33em}2}}}{kC_LT_0}\right).$$
(24)
Esta variação de entropia foi calculada admitindo que a mola está termicamente isolada do ambiente, pelo que é devida a um fenómeno de atrito interno resultante da própria deformação do material.

Se ao contrário admitirmos agora que a mola é um sistema não isolado, de forma a que os valores de temperatura e entropia possam regressar aos seus valores iniciais por troca com o exterior, a energia dissipada na mola acaba por ser transferida integralmente para o ambiente, produzindo neste um aumento $`\mathrm{\Delta }U_0`$=$`W_{diss.}`$. Podemos neste caso estimar o aumento de entropia do ambiente, considerando que este se comporta como uma ”fonte de calor” ($`T_0`$=$`const.`$), usando (12) e (14)

$$\mathrm{\Delta }S_0=\frac{\mathrm{\Delta }U_0}{T_0}=\frac{1}{2}\frac{F_0^{\mathrm{\hspace{0.33em}2}}}{kT_0}.$$
(25)
3 $``$ Entropia e Energia Interna de uma Mola Não Ideal

A análise anterior torna-se um pouco mais complicada no caso de considerarmos uma mola não ideal. Assim, considere-se o exemplo referido em , onde é suposto que num dado intervalo de temperaturas, o comprimento $`L`$ de uma barra, designada aqui em sentido lato por mola, está relacionado com a força $`F`$ que lhe é aplicada através da expressão

$$F=aT^2(LL_0),$$
(26)
onde $`a`$ é uma constante positiva, $`L_0`$ é o comprimento da barra quando não sujeita a qualquer força e $`T`$ é a temperatura da barra quando lhe é aplicada a força $`F`$. Esta barra comporta-se como uma mola não ideal, na medida em que neste caso o termo que correspondia anteriormente à constante de restituição, $`k(T)`$=$`aT^2`$, não é agora constante. Note-se que como iremos mostrar a seguir, mesmo quando a temperatura se mantem constante, e igual à temperatura ambiente $`T_0`$, a mola não se comporta como ideal.

Comecemos por determinar as variações de entropia $`S`$ e energia interna $`U`$ quando a barra é alongada de $`L_0`$ para $`L_1`$ e a temperatura varia entre $`T_0`$ e $`T_1`$. Tal como em , admitamos que a capacidade calorífica para o comprimento $`L_0`$ é dada por $`C_{L_0}`$=$`bT`$, sendo $`b`$ uma constante positiva. Para qualquer outro comprimento tem-se como habitualmente $`C_L`$ dado por (15). Neste caso, o segundo termo de (13), obtido a partir de uma das relações de Maxwell, toma a forma

$$\left(\frac{S}{L}\right)_T=\left(\frac{F}{T}\right)_L=2aT(LL_0).$$
(27)
Calculemos a variação de entropia $`\mathrm{\Delta }`$$`S`$ supondo que a passagem do estado inicial ($`L_0`$,$`T_0`$) para o estado final ($`L_1`$,$`T_1`$), ocorre através de um estado intermédio ($`L_0`$,$`T_1`$). De (13), (16) e (27), podemos escrever

$$(dS)_{L_0}=\frac{C_{L_0}}{T}dT=bdT$$
(28)
e

$$(dS)_{T_1}=2aT_1(LL_0)dL.$$
(29)
Integrando agora (28) e (29), obtém-se para a variação de entropia entre os estados inicial e final

$$\mathrm{\Delta }S=b(T_1T_0)aT_1(L_1L_0)^2.$$
(30)
No que se refere à variação de energia interna $`\mathrm{\Delta }`$$`U`$, ela pode ser calculada diferenciando (30)

$$dS=bdTa(LL_0)^2dT2aT(LL_0)dL$$
(31)
e substituindo esta expressão em (14) usando (26)

$$dU=bTdTaT(LL_0)^2dTaT^2(LL_0)dL.$$
(32)
Uma vez integrada esta última expressão, podemos escrever

$$\mathrm{\Delta }U=\frac{1}{2}b(T_1^{\mathrm{\hspace{0.33em}2}}T_0^{\mathrm{\hspace{0.33em}2}})\frac{1}{2}aT_1^{\mathrm{\hspace{0.33em}2}}(L_1L_0)^2.$$
(33)
É interessante analisarmos antes de mais a que nos conduzem estes resultados numa transformação em que a temperatura se mantenha sempre constante. Assim, vamos supor que a elongação da mola se faz progressivamente, por acréscimos infinitesimais da força exterior, suficientemente lentos para que a mola se mantenha à temperatura ambiente $`T_0`$. Trata-se pois de uma situação onde, de acordo com (30) e (33), as variações de entropia e de energia da mola são dadas por

$$\mathrm{\Delta }S=aT_0(L_1L_0)^2;$$
(34)
$$\mathrm{\Delta }U=\frac{1}{2}aT_0^{\mathrm{\hspace{0.33em}2}}(L_1L_0)^2.$$
(35)
No caso concreto que estamos aqui a tratar, a transformação é reversível, dado que ocorre através de uma sucessão de estados de equilíbrio. Nestas condições, a transformação é isentrópica para o sistema global, constituído pela mola e o ambiente que a rodeia ($`\mathrm{\Delta }S`$+$`\mathrm{\Delta }S_0`$=0), podendo portanto escrever-se de (34)

$$\mathrm{\Delta }S_0=\mathrm{\Delta }S=aT_0(L_1L_0)^2.$$
(36)
Por outro lado, o trabalho realizado pela força aplicada (26) na elongação da mola, neste caso em que a transformação é reversível, com $`T`$=$`T_0`$ ao longo da transformação, é dado por

$$\tau =_{L_0}^{L_1}aT_0^{\mathrm{\hspace{0.33em}2}}(LL_0)𝑑L=\frac{1}{2}aT_0^{\mathrm{\hspace{0.33em}2}}(L_1L_0)^2.$$
(37)
Nestas condições a variação de energia interna $`\mathrm{\Delta }U`$ é dada por

$$\mathrm{\Delta }U=\mathrm{\Delta }U_0+\tau ,$$
(38)
sendo $`\mathrm{\Delta }U_0`$ a energia transferida para o ambiente, a qual é neste caso igual a duas vezes o trabalho realizado pela força exterior na elongação da mola

$$\mathrm{\Delta }U_0=2\tau =aT_0^{\mathrm{\hspace{0.33em}2}}(L_1L_0)^2,$$
(39)
verificando-se ainda $`\mathrm{\Delta }U_0`$=$`T_0`$ $`\mathrm{\Delta }S_0`$.

Repare-se que o resultado a que acabámos de chegar é assaz curioso. Não só o trabalho realizado na elongação da mola é integralmente transferido para o ambiente, como também a energia interna é reduzida exactamente da mesma quantidade, a qual é também transferida para o ambiente. Se procurassemos agora analogias com sistemas eléctricos, encontrariamos uma situação em tudo semelhante, quando se afasta entre si duas armaduras de um condensador plano a potenciais constantes. Neste caso, o trabalho realizado pela força exterior necessário ao afastamento das armaduras é integralmente enviado para a bateria, assim como a própria energia electrostática do sistema diminui do mesmo valor sendo essa energia transferida também para a bateria. A bateria comporta-se aqui como o ambiente no caso da mola, recebendo tanto o trabalho realizado pela força exterior, como a diminuição da energia electrostática do condensador. Obviamente o resultado a que se chegou depende da lei de variação com a temperatura (26). Se outra dependência existisse já esta analogia não se verificaria.

O caso que acabamos de estudar corresponde à situação limite de uma transformação infinitamente lenta e onde a mola se comporta como um sistema não isolado. Considere-se agora o caso em que a mola se encontra isolada termicamente do ambiente e onde a sua elongação ocorre de uma forma isentrópica, $`\mathrm{\Delta }`$$`S`$=0. De (30), podemos escrever

$$b(T_1T_0)=aT_1(L_1L_0)^2$$
(40)
e portanto

$$T_1=\frac{bT_0}{ba(L_1L_0)^2}.$$
(41)
Existe agora um aumento da temperatura com a elongação da mola; $`T_1`$$`>`$$`T_0`$ para $`L_1`$$``$$`L_0`$. Substituindo (41) em (33), obtém-se para a variação de energia interna da mola,

$$\mathrm{\Delta }U=\frac{1}{2}a\frac{(L_1L_0)^2}{ba(L_1L_0)^2}bT_0^{\mathrm{\hspace{0.33em}2}},$$
(42)
e tendo de novo em conta (41),

$$\mathrm{\Delta }U=\frac{1}{2}aT_1T_0(L_1L_0)^2.$$
(43)
Entre dois pontos em que a temperature não varie apreciavelmente, tem-se

$$\mathrm{\Delta }U=\frac{1}{2}k^{}(L_1L_0)^2,$$
(44)
com $`k^{}`$=$`aT_1T_0`$$``$$`const.`$ Aparentemente encontramos de novo a situação da mola ideal, pois entre dois pontos próximos $`T_1`$$``$$`T_0`$ e $`k^{}`$ toma um valor próximo da constante de restituição da mola $`k`$=$`a`$$`T^2`$. Contudo, o primeiro termo de (21) é agora desprezável, pelo que a energia interna da mola só possui o termo de energia potencial.

A situação descrita por (43) tem algum interesse, pois permite aparentemente reeencontrar a expressão vulgarmente usada para a energia de uma mola. Contudo, no caso presente, (43) não representa uma energia potencial.

4 $``$ Conclusões

A Mecânica e a Termodinâmica são habitualmente tratadas como dois capítulos praticamente estanques da Física, tornando-se muito difícil a interpretação de certos fenómenos à luz dos dois conceitos. Neste trabalho procurou-se contribuir, usando um sistema extremamente simples, para que esta integração possa vir para a ordem do dia, chamando a atenção para o facto de que quando se tratam sistemas dissipativos, é desejável que uma análise termodinâmica seja feita também em paralelo com o estudo mecânico. Só desta forma podemos entender o que acontece, de facto, a grandezas tão importantes como a energia interna ou a entropia de um sistema.

Ao longo deste trabalho começou-se por mostrar que a energia dissipada, na elongação de uma mola ideal, por aplicação de uma força $`F_0`$ constante, é sempre igual a metade do trabalho necessário para o fazer, independentemente do meio envolvente. Se este for o vácuo, a energia é dissipada internamente na própria mola, produzindo uma elevação de temperatura, a qual permanecerá ou não indefinidamente na mola consoante esta estiver ou não isolada do ambiente. No caso de uma mola ideal, foram calculados os aumentos de entropia da mola e do ambiente, nos casos em que a mola está, respectivamente, isolada e não isolada termicamente do exterior.

Por último, foi considerado ainda o caso de uma mola não ideal, em que a constante de restituição da mola é uma função da temperatura. Aqui duas situações distintas foram analisadas: Sistema não isolado com a elongação da mola a ter lugar através de um processo infinitamente lento; Sistema isolado com a elongação da mola a ocorrer de uma forma isentrópica. Os aspectos mais relevantes para cada situação foram devidamente evidenciados e discutidos.

$`[`$1$`]`$ $``$ Laufer G., Am. J. Phys. 51, (1983) 42.
$`[`$2$`]`$ $``$ Barrow G. M., J. Chem. Educ. 65, (1988) 122.
$`[`$3$`]`$ $``$ Reif F., Fundamentals of Statistical and Thermal Physics (McGraw-Hill, Tokyo, 1965), p.196.

Agradeço ao Prof. Jorge Loureiro do IST todas as críticas, sugestões e revisão cuidadosa deste manuscrito.


text:
# I Introduction

## I Introduction

The study of neutron-proton pairing especially in the T=0 channel is a particularly prominent topic these days. While the number of journal articles are far too numerous to reference one might begin to make some headway into the varied approaches by starting from the references found in Refs. . In so doing one will find a field of study filled with disagreement and occasionally strife.

The approach here is to study the effects of the T=0 portion of the neutron-proton interaction by removing it altogether. In this way it is hoped that we can understand where in nuclear structure this portion of the interaction is found to play a large role. This would in turn suggest where T=0 pairing might most clearly reveal its presence or absence. In this work we will examine the yrast spectra of the even-even $`(fp)^4`$ nuclei 44Ti,46Ti, 48Cr, and $`{}_{}{}^{50}Cr`$ with two sets of interaction parameters. First we have the full FPD6 interaction , and then we use the same interaction for the two-body T=1 matrix elements while setting all the T=0 two-body interaction matrix elements to zero. We shall denote this interaction as T0FPD6. (This modification of an effective interaction is along the same lines of that used by Satula et al. to examine Wigner energies a few years ago .) We have used this modification of FPD6 in the past to study a variety of things and in particular the full fp spectrum of 44Ti . We found to our surprise and to the surprise of many others that one could obtain a fairly decent spectrum for the levels of 44Ti with the second interaction here labeled T0FPD6. This is surprising because the T=0 interaction is by no means small. The only bound two nucleon system is the T=0 combination of two nucleons, the deuteron. Moreover the largest valued matrix elements in the FPD6 interaction tend to be those in the T=0 channel.

With the one nucleus 44Ti shown in a full fp space calculation in Ref.  we have the first clue that the spectra of even-even nuclei are relatively insensitive to the T=0 two-body interaction matrix elements. To make this more conclusive however we need to examine more nuclei. We have expanded the examination of the even-even nuclei to those listed above. The sample we have chosen consists of nuclei with the same collective properties - ground state bands which have some rotational properties but are not extremely rotational.

Also we now look at transition rates - perhaps these will prove more sensitive to the presence of T=0 two-body matrix elements than the spectra. This will be our interest in a few sections.

It should be noted that in Refs.  a wide range of topics is addressed beyond the spectra of even-even nuclei. These topics include a partial dynamical symmetry that arises when one uses the T0FPD6 interaction in a single j shell for 43Sc and 44Ti. Also while using the T0FPD6 interaction a subtle relationship between the T=$`\frac{1}{2}`$ states in 43Sc and T=$`\frac{3}{2}`$ states in 43Ca likewise between the the T=$`0`$ states in 44Ti and T=$`2`$ states in 44Ca. We also considered even-odd nuclei and addressed the topic of how the T=0 two-body matrix elements affect B(M1) transitions - both spin and orbital components, and Gamow-Teller transitions. In many cases the transition rates were very sensitive to the presence or absence of the T=0 matrix elements. This was especially the case for some orbital B(M1)’s and the Gamow-Teller transitions.

Here things will be kept simple and we focus on the spectra and B(E2)’s of the yrast levels of selected even-even nuclei. We will examine the sensitivity of these observables on the T=0 two-body interaction matrix elements by setting them to zero and comparing the results thus obtained with those when the T=0 matrix elements are reintroduced.

## II Results

As mentioned in the introductions we perform calculations of even-even nuclei with and without the T=0 two-body matrix elements of the FPD6 interaction. Thus each of the figures 1 to 6 corresponding to the nuclei 44Ti, 46Ti, 48Cr, and 50Cr will consist of three columns. The first column is the Yrast spectra calculated with the full FPD6 interaction, the second column with the T0FPD6 interaction, and the third column shows the levels from experiment.

We have previously discussed 44Ti so we will start with 46Ti, but before dissecting the details of the spectra, note the exceptional results in column 2 of figure 2. There we see that the results found when all the T=0 two-body interaction matrix elements are set to zero (T0FPD6) in this complete fp calculation that the resulting spectra looks quite reasonable in comparison with both the full FPD6 and the known experimental levels. Indeed it is difficult to choose between T0FPD6 and FPD6 as to which yields a better fit to experiment.

A closer look shows some differences. The odd spin excitation energies come down by about 1 MeV when the T=0 matrix elements are set to zero. Experimentally some odd spin excitation energies are known (J=1+ and J=11+) but not as many as the even spins. The even spin spectrum is slightly more spread out in the full FPD6 calculation that is to say there is more of a tendency towards a rotational spectrum. The low lying spectrum is fit somewhat better with the full FPD6 interaction but there is a better fit to the J=14+ state using the T0FPD6 interaction.

There are similar stories for 48Cr (Figures 3,4) and 50Cr (Figures 5 and 6). (These nuclei are of some interest as they both display backbends .) In 44Ti there is a bigger difference between FPD6 and T0FPD6 perhaps results from the greater deformation in 44Ti.

As is the case for 46Ti and indeed many of the even-even nuclei in the region very few odd spin states are known. If these were to be made available experimentally it might be easier to demonstrate more clearly the preference between the full FPD6 or the modification T0FPD6. The fact that one can even think of offering into competition an interaction in which all the T=0 two-body interaction matrix elements are set to zero is quite remarkable and perhaps a bit disturbing.

## III B(E2) rates

The B(E2) rates for 44Ti, 46Ti, 48Cr, and 50Cr are listed in tables I to IV. We allow up to t nucleons to be excited from the f7/2 shell to the rest of the fp shell. The values of t used are 4,3,2, and 2 respectively. The effective charges used are the standard $`1.5e`$ for the proton and $`0.5e`$ for the neutron. The difference in the effective charges from 1 and 0 is intended to take care of the fact that the $`\mathrm{\Delta }N=2`$ and higher excitations are not present in this model space. The results for FPD6 and T0FPD6 are shown. We also display the ratios of the results for the 2 interactions.

For 46Ti the reintroduction of the T=0 two-body matrix elements causes an increase (relative to T0FPD6) of a factor of two or more for all the transitions considered. So there is evidence here that the T=0 matrix elements contribute to the collectivity. This is not seen by looking at just the excitation energies in 46Ti where FPD6 and T0FPD6 give very similar results.

In 48Cr and 50Cr a more interesting behavior evolves. In t=2 calculations, the B(E2)’s for low spin transitions get and enhancement with the reintroduction of the T=0 two-body matrix elements but for some of the higher spins we get a suppression. The suppression occurs when backbending occurs as in 50Cr in the J=8 $``$ 10 transition.

In summary, in studying the problem of the T=0 neutron-proton interaction in a nucleus, it may prove more fruitful to begin by removing this channel altogether as was done here by setting all the T=0 two-body matrix elements elements to zero and then reintroducing them, rather than adopting the more common approach of investigating the effects of a pairing interaction separated from the rest of the interaction. This may be especially true in the shell model as the suggestion has been made by Satula and Wyss that it may not be appropriate to separate out a pairing interaction from the rest of the Hamiltonian in a shell model context .

This work was supported by the U.S. Dept. of Energy under Grant No. DE-FG02-95ER-40940 and one of us by GK-12 NSF9979491 (SJQR).


text:
# High Resolution X-Ray Spectroscopic Constraints on Cooling Flow Models of Clusters of Galaxies and Gaseous Haloes Around Elliptical Galaxies

## 1. Introduction

The study of cooling flows in clusters of galaxies and the gaseous haloes of ellipticals has been a controversial area of research for nearly 30 years. As first emphasized independently by Cowie and Binney (1977), Fabian and Nulsen (1977), and Matthews and Bregman (1978), the cooling time for gas at the cores of these systems is often significantly less than the Hubble time. This implies that the gas should continually lose energy by radiation. A pressure driven cooling flow is thought to develop, causing the gas in the densest regions to slowly accrete onto the central dominant galaxies. The net result is a steady deposition of cold matter in the core of the cluster or massive galaxy.

General support for this picture came from some of the first imaging and spectroscopic observations obtained with the Einstein Observatory. Images of many clusters indicate sharp peaks in the surface brightness distribution, as predicted by cooling flow models. Where crude spectrophotometry has been available, one typically sees a marked softening of the spectrum in the central-most regions, again consistent with expectations. Spectroscopically, the presence of the cooling flow should manifest itself as an excess of emission (over isothermal models) at the lowest energies, or via the presence of soft X-ray emission lines inconsistent with the background gas temperature. The latter effect was detected with the Focal Plane Crystal Spectrometer on Einstein, for the gaseous halo around M 87 (Canizares et al., 1979, 1982).

Nevertheless, subsequent observations led to refinements of the simplest cooling flow models, and also pointed to some inconsistencies. For example, for a homogeneous cooling flow, the mass deposition rate should be even more centrally peaked than is inferred from the X-ray brightness profiles (Johnstone et al. 1992). Typically, a mass deposition rate roughly proportional to the radius is derived. This was potentially explained by Nulsen (1986) as a consequence of thermal instability. Dense blobs form locally everywhere in the intracluster medium where the cooling time is less than the age of the system. The result is that the cooling gas resides in a multiphase medium, distributed over a rather large volume. However, this picture was criticized by Balbus and Soker (1989) and others, who argued that in the presence of the gravitational potential, the growth rates of the linear instability are too weak to allow the cooling blobs to form.

A potentially more serious problem is that the end products of cooling are generally not seen in longer wavelength bands at the expected levels. For example, only a very small fraction of the cooling gas can form stars. H I absorption measurements have found no evidence for cold condensed clouds (O’Dea et al. 1998), and while CO emission has been detected (Edge 2001), the amount of inferred molecular gas is a factor 10 below what is predicted by X-ray cooling estimates.

Further, the first moderate resolution spectra obtained with the Einstein Solid State Spectrometer and with ASCA did not provide good quantitative fits to the simplest cooling flow models. Although cool gas is clearly present in the spectrum, there is a dearth of emission at the lowest energies. This was interpreted as evidence for cold absorbing material intermixed with the X-ray emitting gas (White et al. 1991). However, if cold gas is indeed present in the intracluster medium, it remains unclear why it does not form stars.

## 2. The Importance of High Resolution Soft X-ray Spectroscopy

Independent of the details of the cooling flow process, the emitted X-ray spectrum of the cooling gas can be robustly predicted using simple thermodynamic arguments: For an isobaric cooling flow, the differential luminosity emitted in a small temperature range $`dT`$ is given by:

$$dL=\frac{5}{2}\frac{\dot{M}}{\mu m_p}kdT.$$

$`dL`$ can also be expressed in terms of the cooling function of the plasma:

$$dL=\mathrm{\Lambda }(T)dEM$$

where $`dEM=n^2dV`$ is the differential emission measure. Thus:

$$\frac{dEM}{dT}=\frac{5}{2}\frac{\dot{M}}{\mu m_p}\frac{k}{\mathrm{\Lambda }(T)}.$$
(1)
The differential emission measure distribution predicted by this simple argument is displayed in Figure 1. In calculating $`\mathrm{\Lambda }(T)`$, we have assumed one-third cosmic abundances, which is typical of intracluster gas. Marked in the figure are the respective temperatures where each of the Fe L-shell ions becomes the dominant charge state. As can be seen, the Fe L complex very nicely samples the expected distribution. Hence by measuring the strengths of the various Fe L emission lines in a high resolution soft X-ray spectrum, we can accurately constrain cooling flow models.

## 3. Relevance of XMM-Newton

The first experiment capable of unambiguously making such measurements is the Reflection Grating Spectrometer (RGS) on the XMM-Newton Observatory. XMM-Newton incorporates three densely nested grazing incidence telescopes, collectively providing unparalleled area for imaging and spectroscopy. The RGS consists of arrays of grazing incidence reflection gratings, mounted immediately behind two of the three telescopes, which pick off roughly half the light in the beams, and disperse it to sets of dedicated charge-coupled device detectors (CCDs). For each of the two arrays, the gratings are mounted on an inverted Rowland circle, which also includes the respective telescope focus and the readout strip of CCDs. The gratings are all identical, and they are mounted at the same graze angle with respect to the incident ray passing through each grating center. This configuration produces nearly stigmatic and aberration-free focussing at all wavelengths in the spectrum. For a more complete description of this experiment, see den Herder et al. (2002). The other experiments on XMM-Newton are the European Photon Imaging Camera (EPIC), arrays of CCDs mounted at the foci of the three telescopes (see Turner et al. 2002 and Strüder et al. 2002), and the Optical Monitor (OM), an optical/ultraviolet telescope co-aligned with with the X-ray instruments (Mason et al. 2002).

The RGS bandpass ($`\lambda =5`$ to $`38\mathrm{\AA }`$), was adjusted so as to optimally sample the broad range of Fe L-shell transitions, as well as the K-shell features associated with C, N, O, Ne, Mg, and Si. These are the most prominent emission lines expected for collisionally ionized gas, emitting in the temperature regime $`T=10^6`$ to $`10^8K`$, characteristic of the expected cooling flows in clusters and elliptical galaxies. The effective area over this band averages $`120\mathrm{cm}^2`$, if we include both of the two identical instruments. For on-axis point sources, the spectral resolution is $`\mathrm{\Delta }\lambda 0.06`$ Å, varying only slowly as a function of wavelength. Like the transmission grating experiments on the Chandra Observatory, the RGS is a slitless spectrometer, which means that the resolution degrades for extended sources. Nevertheless, since the dispersion of the RGS gratings is very high, the degradation is fairly weak:

$$\mathrm{\Delta }\lambda 0.1\mathrm{\AA }\times \mathrm{\Delta }\theta $$

where $`\mathrm{\Delta }\theta `$ is the source size in arc-minutes. Thus, for only moderately extended sources like clusters and ellipticals, we still obtain fairly high resolution spectra. (In contrast, the spectral resolution of the Chandra experiments degrades by almost a factor of 100 for a source extent of 1 arc-minute.)

## 4. RGS Observations of Clusters and Groups

As a component of its guaranteed time observing program, the RGS instrument consortium observed a sample of 14 compact clusters and groups of galaxies, chosen to exploit the spectral sensitivity of the experiment (see Table 1). For a complete discussion of the analysis of these data sets, see Peterson et al. (2002). The sample includes clusters and groups at a range of temperatures, from 1 to 10 keV. All of these systems were expected to host cooling-flows, with inferred mass deposition rates varying by three orders of magnitude ($`\dot{M}1`$ to $`1000M_{}\mathrm{yr}^1`$). Surprisingly, it was apparent from even the first of these observations, that our data were incompatible with the robust spectral predictions of the standard cooling flow model. An example is illustrated in Figure 2, which displays the derived spectra for the massive cluster A 1835. The intracluster medium in this system has a background temperature $`kT_08`$ keV, and an inferred mass deposition rate $`\dot{M}2300M_{}\mathrm{yr}^1`$ (see Peterson et al. 2001 for a more complete description of the analysis for this case.) In the top panel of the figure, we show the simulated cooling flow spectrum that we expected for the parameters appropriate to this source. The middle panel shows the comparison with the data. As can be seen, the highest ionization lines expected (Fe XXIII, Fe XXIV) are indeed detected in the data, however the model vastly overpredicts the intensity of the lower ionization Fe L-shell lines, as well as the O VIII and O VII K-shell transitions. The discrepancy indicates a dearth of cooler gas. In the bottom panel, we show the comparison to the same cooling flow model, but now with the predicted emission measure distribution truncated at $`kT=2.7`$ keV, with no contributions from gas at lower temperatures. This somewhat arbitrary modification of the model does appear to provide a fairly good description of the data.

Note that the data do not indicate the complete absence of a cooling flow. In particular, the detected Fe XXIII and XXIV lines are indicative of temperatures less than 4 keV, and cannot arise in the 8 keV gas that characterizes most of the intracluster medium. Interestingly, the intensities of these lines are roughly at the right level for the mass deposition rate inferred from the X-ray surface brightness. It is as if the cooling does occur as predicted, but the emergent luminosity somehow disappears from the soft X-ray band at temperatures below $``$ half the background temperature of the cluster.

This same pattern is also found in all of the other clusters in our sample. A subset of the measured spectra are shown in Figure 3. Note that in all cases, Fe XXIII and Fe XXIV lines are detected prominently, while emission lines of lower ionization species (e.g. Fe XVII) are strongly suppressed. While O VIII Lyman-$`\alpha `$ ($`\lambda =18.97\AA `$) is always detected, its measured intensity is also always far weaker than predicted by the cooling flow models for reasonable values of the oxygen abundance (O/Fe $`>`$ 0.2).

## 5. Quantitative Spectral Fits

To derive more quantitative constraints on the discrepancies with the cooling flow predictions, we have performed model fits to the data. Since clusters are extended sources, standard model fitting techniques cannot be easily employed with RGS data. Instead, we use a novel Monte Carlo technique developed expressly for this application (Peterson, Jernigan, and Kahn 2002). The Monte Carlo correctly accounts for the off-axis behavior of the instrument response, and the effects of arbitrary selection cuts and transformations of the data.

We generate photons within the cluster assuming a given spatial and spectral model. These are propagated through the instrument, yielding a distribution in the dispersion coordinate, the cross-dispersion coordinate, and the CCD pulse height. After applying the same cuts as are applied to the real data, the extracted spectrum from the Monte Carlo is compared to the data via a $`\chi ^2`$-statistic. An iterative technique is used to optimize the fit through the adjustment of the astrophysical parameters that specify the spatial/spectral model.

We use a simple $`\beta `$-model to characterize the distribution. The parameters are derived from the simultaneous EPIC images of the source. For the spectrum, we use the MEKAL (Mewe, Kaastra, and Liedahl 1995) model to generate the line and continuum emissivities. The shape of the spectrum is then determined by the differential emission measure distribution.

Outside the “break radius” determined from the EPIC image, we assume that the gas is isothermal, and fit for the temperature and emission measure. Inside, we add in a cooling flow differential emission measure distribution, but we allow the mass deposition rate, $`\dot{M}`$, to be normalized separately in different temperature bins: ($`\frac{1}{2}T_0`$ to $`T_0`$), ($`\frac{1}{4}T_0`$ to $`\frac{1}{2}T_0`$), ($`\frac{1}{8}T_0`$ to $`\frac{1}{4}T_0`$), etc., where $`T_0`$ is the background temperature. This choice of binning is arbitrary, but it works fairly well because the fractional ionization curves are roughly equally spaced in the logarithm of the temperature. The elemental abundances are also left as free parameters, but are assumed to be spatially uniform. The latter assumption does not have a big effect on the fits.

In general, the fits to the data with this model are quite good. This can be seen from Figure 3, where we show a representative sample of the data on individual clusters with the best fit models overlayed. The models appear to correctly account for the relative strengths of the Fe L lines, O VIII, Ne X, and Mg XII.

In Figure 4, we plot the ratio of the derived differential luminosity for each of the temperature bins to the differential luminosity predicted by the simple cooling flow model. In Figure 4a, these points are plotted as a function of absolute temperature, and in Figure 4b, as a function of fractional temperature with respect to the background temperature for each individual cluster. Different clusters are represented by different color points, and the relative points for a few particular clusters are connected by lines for clarity. At the lowest temperatures, we typically derive upper limits, as indicated, but in all cases, we do derive clear detections of emission in at least two distinct temperature bins.

The horizontal dashed lines in these figures represent what would be expected if the cooling flow model were correct. While it is possible for the highest temperature points to lie above the line (due to contributions from the background temperature gas), one would expect the detections and/or upper limits to be consistent with this line for all of the lower temperature points. As can be seen, in essentially all cases, at least some of the lower temperature points lie well below the line, indicating a dearth of cooler gas. When plotted as a function of absolute temperature (Figure 4a), there is no other obvious pattern to the data. However, when plotted as a function of fractional temperature (Figure 4b), it is clear that all of the clusters behave in roughly the same way: there is a clear “drop-out” of emission at temperatures below $`\frac{1}{2}T_0`$, where $`T_0`$ is the background temperature.

In fact, the various clusters appear to be consistent with a differential emission measure distribution given by:

$$\frac{dEM}{dT}=\frac{5}{2}\frac{\dot{M}}{\mu m_p}\frac{k}{\mathrm{\Lambda }(T)}\left(\alpha +1\right)\left(\frac{T}{T_0}\right)^\alpha $$
(2)
with $`\alpha `$ 1 to 2, as opposed to $`\alpha =0`$, as required by Equation (1). This empirical relationship is indicated by the dashed sloping lines in Figure 4b. Thus, in addition to the model being inconsistent with the data at the lowest temperatures, the entire shape of the differential emission measure distribution is inconsistent with the standard cooling-flow model.

## 6. Possible Explanations for the Deficit of Soft Emission

Since the initial reports of some of the results we have presented here (Peterson et al. 2001, Kaastra et al. 2001, Tamura et al. 2001), a number of explanations have been proposed to account for the observed deficit of soft X-ray emission. Most of these work energetically, i.e. they incorporate additional heating or cooling mechanisms that can plausibly supply or remove energy at a rate comparable to the missing soft X-ray luminosity. However, it still unclear whether any can successfully account for the systematic aspects of the observed discrepancies, e.g. the fact that the drop-out of emission always seems to occur at $`\frac{1}{2}T_0`$ over a wide range in $`T_0`$. Below, we briefly discuss a few of the ideas which have been suggested:

Heating. Heating by particle outflows associated with embedded active galactic nuclei (Rosner & Tucker 1989, Tabor & Binney 1993, David et al. 2001), or by subcluster mergers (Markevitch 2001), can supply sufficient energy deposition to counteract the expected cooling. In the context of the multiphase model, the energy input must be distributed evenly everywhere over the thermally unstable region. Otherwise local “pockets” of gas may still cool and condense. In addition, to suppress cooling in a temperature regime that is a fixed fraction of the background temperature, requires that the heating process be self-regulating, i.e. the heating rate has to nearly balance the cooling rate for every cluster, over a range of three orders of magnitude in cooling luminosity. This seems rather contrived, and, at the very least, such models require a fair degree of fine tuning to be considered viable.

Non-standard cooling. A more efficient means of cooling than radiation in the soft X-ray regime could reduce the expected line emission. Simply varying elemental abundances does not help much, since iron dominates both the cooling and the observed spectrum. Mixing the hot gas with cooler embedded gas to reradiate the energy at longer wavelengths (Begelman & Fabian 1990, Fabian et al. 2001) is a possibility, but it is by no means clear why this would always kick in at a fixed fraction of the background temperature. In addition, copious soft X-radiation can be produced by charge transfer at hot/cold interfaces, so it might not be that easy to truly suppress the soft X-ray emission line flux in such a scenario.

Electron conduction. Heat conduction from hotter to colder phases could suppress further cooling as the dense condensations first begin to form (Tucker & Rosner 1983, Stewart et al. 1984, Bertschinger & Meiksin 1986). The presence of tangled magnetic fields in the plasma can strongly suppress conduction down to levels where it becomes unimportant energetically (Chandran & Cowley 1998), but this issue continues to be debated theoretically (Narayan & Medvedev 2001). Observationally, electron conduction has been observed to be suppressed by factors $``$ 100 or more at cold fronts in clusters (Markevitch 2000), so the viability of this explanation also remains suspect.

## 7. RGS Spectroscopy of the Elliptical Galaxy NGC 4636

NGC 4636 is an E0-1 giant elliptical located in the skirt region of the Virgo cluster. The gaseous halo associated with this galaxy has a total X-ray luminosity $`L_X5.7\times 10^{41}\mathrm{ergs}\mathrm{s}^1`$, making it one of the three brightest ellipticals in the sky. The X-ray temperature inferred from ASCA measurements is $`kT0.76`$ keV, with inferred elemental abundances $``$ one third solar (Matsumoto et al. 1997). At this temperature, and the implied electron density, the gas at the core of the halo should cool on timescales less than $`10^8`$ years, leading to a cooling flow with a mass deposition rate $`\dot{M}2M_{}\mathrm{yr}^1`$.

The RGS spectrum of NGC 4636 is shown in Figure 5. As indicated, it is dominated by very strong K-shell line emission from N, O, Ne, and Mg, as well a rich Fe L-shell spectrum. In contrast to the cluster spectra, here we see very strong emission from the lower iron charge states, Fe XVII - XIX. This is a consequence of the much lower background temperature. The red line in the figure is our best fit thermal model for this spectrum, which involves a mean temperature $`0.6`$ keV and elemental abundances ranging from 0.5 to 1 in solar units. For further details on the spectral fit, see Xu et al. 2002.

For this case, however, we can derive additional tight constraints on physical conditions in the halo core by examining the spatial profiles of the data in the cross-dispersion direction. In particular, the RGS provides a nearly stigmatic image of the source in the cross direction. The emission line images are observed to be extended, and the particular profile that we observe for each line provides a clean indication of the distribution of the emission in that line within the halo.

This is of interest for NGC 4636, because, for the inferred halo parameters, the medium should be partially optically thick to resonance scattering for the high oscillator strength Fe L transitions. We observe this effect directly in the data. In particular, the ratio of the Fe XVII $`3s2p`$ lines at 17.1 Å to the Fe XVII $`3d2p`$ lines at 15.0 Å is found to be centrally peaked. Both of these sets of lines have similar emissivity dependences on temperature, so this cannot be produced by a temperature gradient. It is due to the fact that the high oscillator strength 15.0 Å line scatters several times before exiting the halo, thereby broadening its spatial distribution.

In fitting the data, we construct a Monte Carlo spatial/spectral model that is self-consistently constrained to fit both the spectrum and the cross-dispersion emission line profiles. The model is highly over-constrained. For example, the same density required to account for the line intensities, also yields the correct Fe L line optical depths to reproduce the observed spatial profiles. This can only work at one value of the distance, and the distance we infer is consistent with the known distance to the galaxy.

So for this particular case, there is essentially no ambiguity about the physical conditions in the emitting plasma. Nevertheless, the observed spectrum is still dramatically inconsistent with the cooling flow predictions. Specifically, at the implied mass deposition rate, the Fe XVII line intensities should be a factor 10 brighter than observed! Similarly, the O VII line intensities are vastly over-predicted.

Evidently, the multiphase cooling flow model also fails on the scale of an isolated elliptical galaxy, and one without a currently active nucleus. In our view, this single observation provides one of the most stringent challenges to our understanding of cooling in low density cosmic plasmas.

## 8. Summary and Conclusions

The first high spectral resolution X-ray observations of cooling flow clusters continue to challenge conventional models for these systems. Contrary to earlier beliefs, the deficit of soft emission is not due to excess embedded absorbing matter in the intracluster gas, but to a true lack of emission lines from the lowest temperature charge states. We find this problem in all the systems we have studied, from the most massive clusters to the gaseous haloes of individual elliptical galaxies.

Empirically, our data appear consistent with the expression given by Equation (2), where the emergent differential luminosity is proportional to the temperature raised to a power $`\alpha 1`$ to $`2`$, rather than being independent of temperature, as expected from simple thermodynamic arguments. The solution to this puzzle may represent some fundamental new physics in the cooling of low density, hot gas. If so, it may have profound implications for models of galaxy and structure formation in the early universe.

### Acknowledgments.

Analysis of XMM-Newton observations at Columbia University is supported by a grant from NASA. SRON Utrecht is supported financially by NWO, the Netherland Organization for Scientific Research.

## References

Balbus, S. A. & N. Soker 1989, ApJ 341, 611.

Begelman, M. & Fabian, A. C. 1990, MNRAS 244, 26.

Berschinger, E. & A. Meksin 1986, ApJ 306, L1.

Canizares, C. R. et al. 1979, ApJ 234, L33.

Canizares, C. R. et al. 1982, ApJ 262, L32.

Chandran, B. D. G. & S. C. Cowley 1998, PRL 80 Num. 14, 3077.

Cowie, L. L. & J. Binney 1977, ApJ 215, 723.

David, L. P. et al. 2001, ApJ 557, 546.

den Herder et al., A & A 365, 7.

Edge, A. C. 2001, MNRAS 328, 762.

Fabian, A. C. & P. E. J. Nulsen 1977, MNRAS 180, 479.

Fabian, A. C. et al. 2001, MNRAS 320, 20.

Johnstone, R. M. et al. 1992, MNRAS 255, 431-440.

Kaastra, J. S. et al. 2001, A & A 365, 99.

Markevitch, M. et al., 2000, ApJ 541, 542.

Mason, K. O. et al. 2001, A & A 365, 36.

Matthews, W. R. & J. N. Bregman 1978, ApJ 224, 308.

Matsumoto, H. et al. 1997, ApJ 482, 133.

Mewe, R., Kaastra, J., Liedahl, D., 1995, Legacy vol. 6, 16.

Narayan, R. & M. V. Medvedev 2001, ApJ 562, L129.

Nulsen, P. E. J. 1986, MNRAS 221, 377.

O’Dea C. P., Payne, H., & Kocevski 1998, ApJ 116, 623.

Peterson, J. R. et al. 2002, in preparation.

Peterson, J. R. et al. 2001, A & A 365, 104.

Peterson, J. R., J. G. Jernigan, & S. M. Kahn, in preparation; also Peterson, Jernigan, & Kahn, 2002, Proc. SPIE 4847.

Rosner, R. & Tucker, W. 1989, ApJ 338, 761.

Stewart et al. 1984, ApJ 278, 536.

Strüder, L. et al. 2001, A&A 365, 18.

Tabor, G. & J. Binney 1993, MNRAS 263, 123.

Tamura, T. et al. 2001, A & A 365, 87.

Tucker, W. H. & R. Rosner 1983, ApJ 267, 547.

Turner, M. J. L. et al. 2001, A & A 365, 27.

White, D. A. et al. 1991, MNRAS 252, 1991.

Xu, H. et al. 2002, astroph-0110013, ApJ in press.