text:
# Birth and long-time stabilization of out-of-equilibrium coherent structures

## 1 Introduction

The Hamiltonian Mean Field model (HMF) has attracted much attention in the recent years as a toy model to study the dynamics of systems with long-range interactions, and its relation to thermodynamics Antoni ; AnteneodoTsallis ; LatoraRapisardaRuffo . The HMF model describes an assembly of $`N`$ fully coupled rotators, whose Hamiltonian is:

$`H`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{p_i^2}{2}}+{\displaystyle \frac{c}{2N}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\mathrm{cos}(\theta _i\theta _j),`$ (1)
where $`\theta _i`$ is the angle of the $`i`$-th planar rotator with a fixed axis. As the interaction only depends on the angles of the rotators, this model can alternatively be viewed as representing particles that move on a circle, whose positions are given by the $`\theta _i`$ and interact via an infinite-range force. If $`c`$ is negative, the interaction among rotators is ferromagnetic, corresponding to an attractive interaction among particles. When $`c`$ is positive, the interaction among rotators is antiferromagnetic and repulsive in the particle interpretation. In this article, we will focus our study on this latter case, and put $`c=1`$.

As first noticed in Ref. Antoni , this model has a very interesting dynamical behavior. In contrast to the statistical mechanics predictions, a bicluster forms at low energy (see Fig. 1), for a special, but wide class of initial conditions: it appears as soon as the initial velocity dispersion of the particles is small, for any initial spatial distribution. This clustering phenomenon and its unexpected dependence on initial conditions were precisely studied in Ref. drh , but remained essentially unexplained. Several details of the dynamics were numerically studied and it has been in particular emphasized that the energy temperature relation is modified: although still linear, the slope in the molecular dynamics simulations is different from the theoretical prediction. No sign in the numerics was found of the decay to the constant density profile predicted by the canonical ensemble. On the contrary, recently performed simulations FFR  with a smaller number of particles have shown a long-time degradation of the bicluster, suggesting its transient non-equilibrium nature. The question of the time asymptotic stability of the bicluster, in the limit of an infinite number of particles, remains however open and will be discussed in the conclusions.

We will consider in this work two different theoretical approaches that lead to an analytical explanation of the clustering process and of the long time stability of the structure. In section 2, in order to explain the short time bicluster formation, we first propose a hydrodynamical description that links Hamiltonian (1) to a forced Burgers equation (these results were shortly presented in a previous brief note BDR ). This striking phenomenon is then described in Lagrangian coordinates, giving rise to intersections of trajectories and formation of high (or even infinite) particle densities. We then solve in section 3 this equation, using extensively the method of characteristics. Beyond the theoretical interest in such structures for Hamiltonian systems, we will show some connections of this problem with other subjects; namely, active transport in hydrodynamics, and the formation of caustics in the flow of the associated Burgers equation.

Section 4 presents the second approach, which is particularly powerful for explaining the long-time evolution of the bicluster. The dynamics involves two well separated time-scales. Using a variational approach, we will show how to construct an effective Hamiltonian where all the fast oscillations are averaged out. Not only this effective dynamics accurately represents the one of the original Hamiltonian but, in addition, the complete statistical thermodynamics can be easily derived in the microcanonical and canonical ensembles. It predicts the existence of the bicluster as an equilibrium state, with the correct density profile and the microcanonical temperature energy relation found in the numerical experiments of Ref. drh . The long lifetime of these out-of-equilibrium states can therefore be interpreted by the fact that they appear as equilibrium states of an effective Hamiltonian representing the long-time motion. Section 5 is devoted to some conclusions and perspectives of future developments.

## 2 Hydrodynamical description at zero temperature

### 2.1 Introduction

To describe the initial clustering process, we consider the associated Vlasov equation, which can be rigorously derived in the thermodynamic limit $`N\mathrm{}`$ Kinlimit ; BraunHepp . Denoting by $`f(\theta ,p,t)`$ the one particle distribution function, we have:

$$\frac{f}{t}+p\frac{f}{\theta }\left[\frac{1}{2\pi }_{\mathrm{}}^+\mathrm{}𝑑u_0^{2\pi }𝑑\alpha f(\alpha ,u,t)\mathrm{sin}(\theta \alpha )\right]\frac{f}{p}=0$$
(2)
Let us define as follows, a density

$$\rho (\theta ,t)=_{\mathrm{}}^+\mathrm{}f(\theta ,p,t)𝑑p$$
(3)
and a velocity field

$$\rho (\theta ,t)v(\theta ,t)=_{\mathrm{}}^+\mathrm{}pf(\theta ,p,t)𝑑p.$$
(4)
As the numerical simulations reported in Ref. drh  have shown that the bicluster appears when the velocity dispersion is small, we will consider here the zero temperature approximation, i.e. we neglect the velocity dispersion and consider the ansatz $`f(\theta ,p,t)=\rho (\theta ,t)\delta (pv(\theta ,t))`$, where $`\delta `$ refers to the Dirac function. The Vlasov equation (2) can consequently be simplified to

$$\frac{f}{t}+p\frac{f}{\theta }F(\theta ,t)\frac{f}{p}=0$$
(5)
where one notes

$$F(\theta ,t)=\frac{1}{2\pi }_0^{2\pi }𝑑\alpha \rho (\alpha ,t)\mathrm{sin}(\theta \alpha ).$$
(6)
The sum of the time derivative of Eq. (3) and the $`\theta `$-derivative of Eq. (4), leads to

$`{\displaystyle \frac{\rho }{t}}+{\displaystyle \frac{(\rho v)}{\theta }}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\left({\displaystyle \frac{f}{t}}+p{\displaystyle \frac{f}{\theta }}\right)𝑑p`$ (7)
$`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑pF(\theta ,t){\displaystyle \frac{f}{p}}=0.`$
We obtain therefore the equation which accounts for mass conservation

$$\frac{\rho }{t}+\frac{(\rho v)}{\theta }=0.$$
(8)
Using the time derivative of Eq. (4) and Eq. (8), we obtain

$$\frac{\left(\rho v\right)}{t}=\rho \left(\frac{v}{t}+v\frac{v}{\theta }\right)\frac{(\rho v^2)}{\theta }.$$
(9)
However, the time derivative of the left-hand-side gives also

$`{\displaystyle \frac{\left(\rho v\right)}{t}}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}p{\displaystyle \frac{f}{t}}𝑑p`$ (10)
$`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}p^2{\displaystyle \frac{f}{\theta }}𝑑p+F(\theta ,t){\displaystyle _{\mathrm{}}^+\mathrm{}}p{\displaystyle \frac{f}{p}}𝑑p`$
$`=`$ $`v^2{\displaystyle \frac{\rho }{\theta }}\rho {\displaystyle \frac{v}{\theta }}2v+F(\theta ,t)\rho (\theta ,t).`$
Equations (9) and Eq. (10) lead finally to the following Euler equation without pressure term:

$`{\displaystyle \frac{v}{t}}+v{\displaystyle \frac{v}{\theta }}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑\alpha \rho (\alpha ,t)\mathrm{sin}(\theta \alpha )`$ (11)
Equations (8) and (11) show that the dynamics of the model at low temperature can be therefore mapped onto an active scalar advection problem.

### 2.2 Linear analysis

Linearizing equations (8) and (11), by assuming small velocities ($`v1`$) and almost uniform density ($`\rho =1+\rho _1(\theta ,t)`$ with $`\rho _11`$), we are left with:

$`{\displaystyle \frac{\rho _1}{t}}+{\displaystyle \frac{v}{\theta }}=0`$ (12)
$`{\displaystyle \frac{v}{t}}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}\rho _1(\alpha ,t)\mathrm{sin}(\theta \alpha )`$ (13)
a system which is easily solved by a spatial Fourier series development. Assuming a zero initial velocity field, we get

$`\rho _1(\theta ,t)=\sqrt{2}v_m\mathrm{cos}\theta \mathrm{cos}\omega t`$ (14)
$`v(\theta ,t)=v_m\mathrm{sin}\theta \mathrm{sin}\omega t,`$ (15)
where the time-scale of the oscillation is given by the inverse of a sort of “plasma frequency” $`\omega =\sqrt{2}/2`$, which describes, as in the Poisson-Boltzmann case, the instability of the uniform density state Nicholson . In the attractive case ($`c=1`$ in Hamiltonian (1)), this time scale would control the depart from the initial uniform density towards the formation of the density profile of the single cluster that appears at low temperature Ina1 ; Ina2 , in analogy with Jeans instability in gravitational systems Binney .

This linear analysis is, however, not sufficient to explain the formation of the bicluster and we have to carry out a non linear analysis.

### 2.3 Non linear analysis

This analysis relies on the existence of two time-scales in the system: the first one is intrinsic and corresponds to the inverse of the “plasma frequency” $`\omega `$; the second one is connected to the energy per particle $`e`$. When the energy is sufficiently small, these time-scales are very different, and it becomes possible to use averaging methods. The presence of two well separated time-scales also explains, in some sense, the appearance of the bicluster in the low energy limit.

We introduce therefore a long time-scale $`\tau =\epsilon t`$, with $`\epsilon =v_m/\sqrt{2}=\sqrt{2e}`$, and we look for solutions of the following form:

$`v(\theta ,t)`$ $`=`$ $`v_m\mathrm{sin}\theta \mathrm{sin}\omega t+\epsilon u(\theta ,\tau ).`$ (16)
In order to evaluate the force on the r.h.s. of Eq. (11) in the nonlinear regime, we will use expressions (14) and (15), given by the linear analysis. Indeed, due to the special form of the Hamiltonian, the force depends only on the first Fourier component of the density. Hence, our hypothesis amounts to assume that the sinusoïdal behavior of the density, found in the linear regime, holds the same in the non-linear regime: the results presented below confirm the validity of this assumption.

In analogy with studies of the wave-particle interaction in plasma physics plasma ; firpo , our system may be seen as a bulk of particles interacting with waves that are sustained by the bulk itself. Here, the wave is created by the small density and velocity oscillations already present in the linear analysis. In the low energy regime, the phase velocity of the wave with the plasma frequency $`\omega `$ is much higher than the velocities of the bulk particles, causing a very small wave-particle interaction strength: the wave has therefore a very long life-time. This also explains the quality of the linear approximation for the force that we have used in the Euler equation (11).

Introducing expression (16) in Eq. (11), terms to first order in $`\epsilon `$ disappear by construction. Order $`\epsilon ^2`$ terms give:

$`{\displaystyle \frac{u}{\tau }}+u{\displaystyle \frac{u}{\theta }}`$ $`+`$ $`2\mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}^2\omega t`$ (17)
$`+`$ $`\left(\mathrm{sin}\theta {\displaystyle \frac{u}{\theta }}+u\mathrm{cos}\theta \right)\sqrt{2}\mathrm{sin}\omega t=0`$
Averaging over the short time scale $`t`$, we obtain:

$`{\displaystyle \frac{u}{\tau }}+u{\displaystyle \frac{u}{\theta }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}2\theta `$ (18)
which is a spatially forced Burgers equation without viscosity, describing the motion of fluid particles in the potential $`V(\theta )=1/4\mathrm{cos}2\theta `$.

The Burgers equation has been studied in many different contexts by mathematicians interested in a pressureless description of fluid motion, and by physicists in the context of ballistic aggregation models young . In particular, the forced Burgers equation has been very carefully studied from a mathematical viewpoint moserjauslin . Cosmologists have proposed a very interesting adhesion model that uses Burgers-like equations to explain structure formation in the universe youngIHP . Moreover, this equation may be also viewed as a nice model of active transport of the vorticity field in fluid mechanics Castillo .

A well known property of the Burgers equation without viscosity, is that the solution becomes multi-stream after a finite time: the appearance of shocks (see Fig. 2) in the velocity profile $`u(\theta )`$ corresponds to the creation of singularities in the density profile (see Fig. 3). This is a consequence of the two main assumptions: the medium is supposed to be continuous and the temperature is set at zero. Weakening either assumption would have smoothed out the singularities. In the original discrete Hamiltonian model, particles can cross and, after some time, those that travel faster can be catched by the slower ones that lie downstream, creating the spiral dynamics exemplified in Fig. (4). The particle density has a peak at the center of the spiral (Fig. 3), and diverges at the peak position in the $`N\mathrm{}`$ limit at the shock time $`\tau _s`$. The presence of a shock at finite time in the forced Burgers equation does not prevent the hydrodynamical description from remaining valid at longer times, when, moreover, several repeated shocks are observed at regular time intervals. We will indeed show that it can adequately describe even the long-time evolution. In addition, let us remark that it’s the double well shape of the potential that forces the Burgers equation, which is responsible for the formation of two clusters. Starting from the initially uniform state, the particles will move around the bottom of the both wells and the majority of them happen to be there at the same time at the shock time, creating therefore two coherent structures that after several oscillations form the bicluster at long-time.

## 3 Solution of the spatially forced Burgers equation by the method of characteristics

The method of characteristics is the appropriate mathematical tool to study more quantitatively the forced Burgers equation (18), but it has also a direct physical meaning: the characteristics correspond to the Lagrangian trajectories of the Euler equation, and are therefore good approximations for the particle trajectories of the finite $`N`$ Hamiltonian system (1) when $`N1`$.

In the case of unidirectional nonlinear wave motions, the method is standard and proceeds as indicated, for example, in Ref. Zauderer . We obtain in this case the following pendulum equation for the characteristics $`\theta (\tau )`$:

$`{\displaystyle \frac{d^2\theta }{d\tau ^2}}+{\displaystyle \frac{1}{2}}\mathrm{sin}2\theta `$ $`=`$ $`0,`$ (19)
with the initial conditions $`\theta (0)=\theta _0`$ and $`\dot{\theta }(0)=0`$. The solution of the above equation can be easily derived lawden  and reads:

$$\theta =\mathrm{arcsin}\left[\mathrm{sin}\theta _0\text{sn}(\tau +K_0,\mathrm{sin}\theta _0)\right]$$
(20)
where $`K_0=K(\mathrm{sin}\theta _0)`$, $`K`$ is the complete elliptic integral of the 1st kind and sn, the elliptic sine function. Fig. (5a) shows the characteristics, i.e. the solutions of Eq. (19) corresponding to different initial conditions $`\theta _0`$ evenly distributed in $`[0,2\pi ]`$ at $`t=0`$ and with zero initial velocity. One clearly sees their oscillatory behavior inside one of the two wells of the effective potential $`\mathrm{cos}2\theta `$, located at $`\theta =0`$, with a period $`4K_0`$ which strongly depends on the initial position. Fig. (5a) emphasizes also that the characteristics cross themselves, and the associated caustics are responsible for the enhancement of the particle density field that forms the “chevrons”. The time of the first divergence in the density $`\tau _s`$ is the time of the first characteristic crossing and also the time of the first shock in the velocity profile. The particles performing a quasi-harmonic oscillatory motion in the bottom of the well are the ones with shorter periods. The shortest period is obtained in the harmonic limit $`4K(0)=2\pi `$, leading to $`\tau _s=\pi /2`$ for the time of appearance of the first shock. In the original units, it reads

$$t_s=\frac{\tau _s}{\epsilon }=\frac{\pi }{\sqrt{8e}}.$$
(21)
This shock results in a singularity in Eulerian space at the corresponding time. However this singularity disappears immediately after forming and two singularities of another kind arise in its place, which are the boundaries of the three streams region. The chevron is their manifestation in the density profile. The shock time singularity is called $`A_3`$ in Arnold’s classification arnold ; Zeldovich  while the “chevrons” singularity is of $`A_2`$ type. The latter is the boundary of structures and can exist at any moment of time, whereas the former exist only when a chevron originates. The recurrence time for the appearance of chevrons is $`t_n=t_s(2n1)`$ and the comparison with the numerics shows a very good agreement.

Figs. (5a) and (5b) show that the agreement is really good not only at the qualitative but also at the quantitative level. Only two features are missed by the Lagrangian description. The fast oscillations of very small amplitude, already hardly visible in Fig. (5b) is of course totally absent in the characteristics because of the averaging procedure used to get Eq. (18). One can also easily show that the amplitude of this fast oscillation is rapidly decreasing with the increase of the number of particles and is vanishingly small for $`N\mathrm{}`$. The second point is the presence of untrapped particles in Fig. (5b), which is a direct consequence of the oscillations of the height of the potential barriers. This effect is important only for highly energetic particles whose trajectories in phase space are closed to the separatrix of the Kelvin’s cat eye of Fig. (4). However, as we will show in the next section, these particles do not take part to the creation of caustics (i.e. of the chevrons of Fig. 1) which are generated only by the particles close to the bottom of the wells.

### 3.1 The chevrons as caustics of the characteristics

We will now explain the shape of the ”chevrons” (Fig. 1), which correspond to zones of infinite density i.e. to the envelops of the characteristics, the so-called caustics. The family of the characteristics is defined by

$$F(\tau ,\theta ,\theta _0)=\mathrm{sin}\theta (\tau ,\theta _0)\mathrm{sin}(\theta )=0$$
(22)
in the plane $`(\theta ,\tau )`$ with the parameter $`\theta _0`$. The envelope of this family would then be defined by the two following equations

$`F(\tau ,\theta ,\theta _0)`$ $`=`$ $`0`$ (23)
$`{\displaystyle \frac{F}{\theta _0}}(\tau ,\theta ,\theta _0)`$ $`=`$ $`0`$ (24)
However, it is not possible to extract $`\theta _0`$ as a function of $`\theta `$ and $`\tau `$ from Eq. (24), in order to obtain a closed expression $`F(\tau ,\theta ,\theta _0(\tau ,\theta ))`$ for the caustics. We will use approximate expressions of $`F`$ close to the shocks.

In the neighborhood of the shock, the trajectories can be approximated by straight lines, with the following expression:

$`\theta (\tau ,\theta _0)`$ $`=`$ $`v_0\left(\tau K_0\right)`$ (25)
where $`v_0=v(K_0)`$ is the speed of the particle at the bottom of the potential well.

The envelope of this family of characteristics can then be obtained Zauderer  by solving the following system of equations

$`\theta (\tau ,\theta _0)`$ $`=`$ $`v_0\left(\tau K_0\right)`$ (26)
$`0`$ $`=`$ $`v_0^{^{}}\left(\tau K_0\right)v_0K_0^{^{}}`$ (27)
where the primes denotes the derivative with respect to $`\theta _0`$. Using the conservation of energy, we easily get $`v_0=\mathrm{sin}\theta _0`$ and we obtain thus

$`\tau `$ $`=`$ $`\mathrm{tan}\theta _0K_0^{^{}}+K_0`$ (28)
$`\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}^2\theta _0}{\mathrm{cos}\theta _0}}K_0^{^{}}.`$ (29)
It is however possible to go further by using higher order terms in the expression of $`\theta `$ given in Eq. (25), i.e. by considering curves rather than simply straight lines. We have thus

$`\theta (\tau ,\theta _0)`$ $`=`$ $`\theta (K_0)+{\displaystyle \frac{d\theta }{d\tau }}_{|K_0}\left(\tau K_0\right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2\theta }{d\tau ^2}}_{|K_0}\left(\tau K_0\right)^2`$ (30)
$`+{\displaystyle \frac{1}{6}}{\displaystyle \frac{d^3\theta }{d\tau ^3}}_{|K_0}\left(\tau K_0\right)^3+\mathrm{}`$
$`=`$ $`v_0\left(\tau K_0\right){\displaystyle \frac{v_0}{6}}\left(\tau K_0\right)^3+\mathrm{}`$ (31)
where we have replaced in Eq. (31) the derivatives by their expressions. Let us note in particular, that all derivatives of even order will vanish because of the parity property of $`\theta `$. Using the development of $`K_0`$ at the same order, we finally end with the next order approximation for the caustics

$`\tau `$ $`=`$ $`\tau _s+3a\theta _0^2+\left(5b+{\displaystyle \frac{2a}{3}}\right)\theta _0^4`$ (32)
$`+\left(7c+{\displaystyle \frac{4b}{3}}+{\displaystyle \frac{4a}{15}}{\displaystyle \frac{8a^3}{3}}\right)\theta _0^6+O(\theta _0^8)`$
$`\theta `$ $`=`$ $`2a\theta _0^3+\left(4b+{\displaystyle \frac{a}{3}}\right)\theta _0^5`$ (33)
$`+\left(6c+{\displaystyle \frac{2b}{3}}+{\displaystyle \frac{31a}{180}}4a^3\right)\theta _0^7+O(\theta _0^9)`$
where $`(a,b,c)=\frac{\pi }{8}(1,\frac{11}{48},\frac{173}{2880})`$.

It is possible to see in Fig. (6) that the agreement is excellent and that this procedure is really accurate to describe the first chevron. Moreover, one can check that because of the non isochronism of the oscillations, the more energetic particles arrive too late in the bottom of the well to take part to the creation of the caustics, i.e. the singularity in the density space. Consequently, the untrapped particles of the real dynamics, visible in Fig. (5b), but absent of the averaged dynamics of the associated Burgers equation (see Fig. 5a), are irrelevant and do not affect the dynamics of the caustic formation.

This is confirmed by Fig. (7) where the caustics determined by the Burgers’ approach are superimposed on the real trajectories: the shape of the chevrons is very well described. However, the time occurrence of the shock was reduced by a factor 0.84. The reason, clarified in the next section, is due to the renormalization of the oscillation frequency $`\omega `$ which modifies the oscillation period of the particles. As shown below using a self-consistent procedure, $`\omega `$ is reduced with respect to the plasma frequency $`\omega _p`$. As the potential height is inversely proportional to $`\omega `$ this leads therefore to a faster appearance of the shocks.

As already mentioned, the time occurrence of the $`n`$th shock is $`\tau _n(\theta _0)=(2n1)\tau _s`$ which means that its shape will be obtained simply by replacing in Eq. (26) $`K_0`$ by $`(2n1)K_0`$. Therefore, the lowest order approximation of the $`n`$th shocks is

$`\theta `$ $`=`$ $`2a(2n1)\left({\displaystyle \frac{\tau \tau _n}{3a(2n1)}}\right)^{3/2}`$ (34)
$`\theta `$ $``$ $`{\displaystyle \frac{\left(t(2n1)t_s\right)^{3/2}}{\sqrt{2n1}}}`$ (35)
where $`n`$ is the number of appearance of the ”chevron”. The $`1/\sqrt{n}`$ factor accounts for the shrinking of the “chevrons” and the bolded line in Fig. (6) attests that this expression is particularly accurate.

Similar caustics are encountered in astrophysics, to explain the large scale structure of the universe: clusters and super clusters of galaxies are believed to be reminiscent of three dimensional caustics arising from the evolution of an initially slightly inhomogeneous plasma Vergassola ; Zeldovich .

### 3.2 Probability distribution

At this point, it would be important to calculate the density distribution in the vicinity of the singularity. For this hydrodynamical approximation, it would be fully determined by the initial conditions, i.e. the initial velocity distribution. Let us recall that we earlier found drh , that the following analytical formula

$$𝒫(\theta )=\frac{1}{2\pi }\left(1\mathrm{log}(2|\mathrm{sin}\theta |)\right)$$
(36)
was an excellent approximation (see Fig. 8), even if we didn’t have any theory for its derivation. Here using the Burgers’ approach, let us obtain a similar result.

Along the particle trajectory, one has

$$\dot{\theta }=\sqrt{\frac{\mathrm{cos}2\theta \mathrm{cos}2\theta _0}{2}}$$
(37)
The time $`dt(\theta )`$ spent by a particle close to a position $`\theta `$ is inversely proportional to its velocity $`\dot{\theta }(\theta _0)`$. Therefore the trajectory parametrized by $`\theta _0`$ (and initiated in $`\theta _0`$) gives a contribution to the density

$`\rho _{\theta _0}(\theta )d\theta `$ $`=`$ $`{\displaystyle \frac{dt(\theta )d\theta }{{\displaystyle _{\theta _0}^{\theta _0}}𝑑t(\theta )}}={\displaystyle \frac{2\sqrt{2}d\theta }{T(\theta _0)\sqrt{\mathrm{cos}2\theta \mathrm{cos}2\theta _0}}}`$ (38)
where $`T(\theta _0)`$ is the period of the trajectory. By recalling that the initial distribution is homogeneous in these numerical computations, one obtains the probability density by averaging over the time, the density of characteristics at a fixed position. We obtain:

$`\rho (\theta ){\displaystyle _{2\theta }^\pi }{\displaystyle \frac{2\sqrt{2}du}{T(u)\sqrt{\mathrm{cos}2\theta \mathrm{cos}2u}}}`$ (39)
for $`\theta =[0,\pi /2]`$ whereas the whole distribution is obtained by symmetry and $`\pi `$-periodicity.

Due to the numerous approximations, the agreement for long time is not perfect, but equation (39) gives a good result as shown by Fig. (8). The disagreement is visible in $`0`$ and $`\pi `$, i.e. close to the separatrix. Actually, as we will see in the following, the amplitude of the effective potential slightly oscillates; this allows the existence of untrapped particles which smooth the density.

As a conclusion, this hydrodynamical approach has proved to be very successful, since it explains qualitatively and quantitatively, for short times, the main features of the bicluster formation. Nevertheless, the stability and the long time behavior of the structure is not yet understood: this is the object of the next section.

## 4 Long time evolution of the bicluster

Since the thermodynamics of the antiferromagnetic HMF model predicts a perfectly homogeneous equilibrium state, the bicluster is a non equilibrium structure. However, its long time stability, as well as the fact that it is reached starting from a variety of initial conditions suggest some statistical explanation.

In addition, we have learned from the previous section that particle motion may be split into two parts: a small amplitude rapid oscillation superimposed to a large amplitude slow motion. The idea is thus to look for an effective Hamiltonian dynamics which would only retain the slow motion, and would be appropriate for a statistical description.

### 4.1 The fast variables

The complete Lagrangian of the antiferromagnetic HMF being

$`L(\theta _i,\dot{\theta _i})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\dot{\theta _i}^2}{2}}{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i,j}{}}\mathrm{cos}(\theta _i\theta _j),`$ (40)
the equations of motion are obtained by minimizing the action $`S=L𝑑t`$ with respect to the functions $`\theta _i`$. Taking advantage of the two well defined time scales, we introduce the following ansatz:

$$\theta _i=\theta _i^0(\tau )+\epsilon f_i(t,\tau )$$
(41)
where $`\tau =\epsilon t`$ and $`\epsilon =\sqrt{2e}`$ as already defined in the section (2.3). The $`f_i`$’s represent the small rapid oscillations and the $`\theta _i^0`$ the slow variables which we are finally interested in. The idea is then to use a variational approach (least action principle), first to obtain the equations of motion of the fast variables $`f_i`$, and then to reintroduce the solutions into the action. In a second stage, averaging over the fast time $`t`$, we will obtain a variational principle for the slow variables. This approach, inspired by witham , has the double advantage of allowing to obtain an Hamiltonian system for the slow variables, and the one of exhibiting a natural conservation law for this dynamics (the adiabatic invariant). Appendix A briefly presents the use of method witham  for a slowly modulated harmonic oscillator, in order to emphasize its power for a very simple model, away from the rather technical context presented here.

We first notice that since the energy of the system is by definition of order $`\epsilon ^2`$, the sums $`1/N\mathrm{cos}\theta _i^0`$ and $`1/N\mathrm{sin}\theta _i^0`$, which represent the two components of the magnetization vector in terms of the slow variables $`\theta _i^0`$, are of order $`\epsilon `$. We thus define

$`\epsilon M_{1x}^0`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\mathrm{cos}\left(\theta _i^0+\psi \right)`$ (42)
$`\epsilon M_{1y}^0`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\mathrm{sin}\left(\theta _i^0+\psi \right)`$ (43)
where the phase $`\psi `$ $`[0,\pi ]`$ is chosen such that the scalar dynamical indicator of the clustering is

$$|M_2|=\left|\frac{1}{N}\underset{i}{}\mathrm{exp}\left(2i\theta _i^0\right)\right|=\frac{1}{N}\underset{i}{}\mathrm{cos}2\left(\theta _i^0+\psi \right).$$
(44)
We now introduce the ansatz (41) into the Lagrangian (40), and develop the cosine up to order $`\epsilon ^2`$, obtaining the new Lagrangian $`L_2`$ that depends on the $`\theta _i^0`$’s, the $`f_i`$’s and their time derivatives

$`L_2`$ $`=`$ $`L_2(\theta _i^0,{\displaystyle \frac{d\theta _i^0}{d\tau }},f_i,{\displaystyle \frac{df_i}{dt}})`$ (45)
$`=`$ $`{\displaystyle \frac{\epsilon ^2}{2}}{\displaystyle \underset{i}{}}\left[\left({\displaystyle \frac{d\theta _i^0}{d\tau }}\right)^2+2{\displaystyle \frac{d\theta _i^0}{d\tau }}{\displaystyle \frac{df_i}{dt}}+\left({\displaystyle \frac{df_i}{dt}}\right)^2\right]`$
$``$ $`{\displaystyle \frac{\epsilon ^2N}{2}}\left[(M_{1x}^0)^2+(M_{1y}^0)^2\right]+{\displaystyle \frac{\epsilon }{N}}{\displaystyle \underset{i,j}{}}f_i\mathrm{sin}(\theta _i^0\theta _j^0)`$
$``$ $`{\displaystyle \frac{\epsilon ^2}{2N}}{\displaystyle \underset{i,j}{}}f_if_j\mathrm{cos}(\theta _i^0\theta _j^0).`$
From this expression, considering $`\tau `$ as a constant, we write the equations of motion for the fast variables $`f_i`$’s. We get

$`\ddot{f_i}`$ $`=M_{1x}^0\mathrm{sin}\left(\theta _i^0+\psi \right)`$ $`M_{1y}^0\mathrm{cos}\left(\theta _i^0+\psi \right)`$ (46)
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j}{}}f_j\mathrm{cos}\left(\theta _i^0\theta _j^0\right).`$
This is a linear second order equation, with constant coefficients with respect to the fast time, whose solution thus requires the diagonalization of the $`N`$x$`N`$ matrix

$$𝐓=\left[T_{ij}\right]=\frac{1}{N}\left[\mathrm{cos}(\theta _i^0\theta _j^0)\right],$$
(47)
which has only two non zero eigenvalues (see Appendix B)

$$\omega _\pm ^2=\frac{1\pm |M_2|}{2}.$$
(48)
The eigenvectors corresponding to $`\omega _\pm ^2`$ are respectively

$`𝐗_+`$ $`=`$ $`\left[\mathrm{cos}(\theta _i^0+\psi )\right]_{i=1\mathrm{}N}`$ (49)
$`𝐗_{}`$ $`=`$ $`\left[\mathrm{sin}(\theta _i^0+\psi )\right]_{i=1\mathrm{}N}.`$ (50)
The general solution for the $`f_i`$’s is therefore

$`f_i(t,\tau )`$ $`=`$ $`\left[\sqrt{2}A_+\mathrm{sin}(\omega _+t+\phi _+){\displaystyle \frac{M_{1y}^0}{\omega _+^2}}\right]\mathrm{cos}(\theta _i^0+\psi )`$ (51)
$`+`$ $`\left[\sqrt{2}A_{}\mathrm{sin}(\omega _{}t+\phi _{})+{\displaystyle \frac{M_{1x}^0}{\omega _{}^2}}\right]\mathrm{sin}(\theta _i^0+\psi ).`$
### 4.2 The slow variables

The previous solution for the $`f_i`$ suggests the following new ansatz

$`f_i`$ $`=`$ $`\left[\sqrt{2}A_+(\tau )\mathrm{sin}(\varphi _+(t))+a_+(\tau )\right]\mathrm{cos}(\theta _i^0+\psi )`$ (52)
$`+`$ $`\left[\sqrt{2}A_{}(\tau )\mathrm{sin}(\varphi _{}(t))+a_{}(\tau )\right]\mathrm{sin}(\theta _i^0+\psi )`$
where $`\varphi _\pm (t)`$ are fast variables, and $`d\varphi _\pm /dt`$, $`A_\pm (\tau )`$, and $`a_\pm (\tau )`$ are slow variables. We introduce (52) in the Lagrangian (45), and we average over the fast variables $`\varphi _\pm `$. Dropping the $`\epsilon ^2`$ overall factor, the averaged Lagrangian reads

$`L_{eff}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{d\theta _i^0}{d\tau }}\right)^2+{\displaystyle \frac{N}{2}}\left[A_+^2\omega _+^2\dot{\varphi }_+^2+A_{}^2\omega _{}^2\dot{\varphi }_{}^2\right]`$ (53)
$`N[{\displaystyle \frac{M_{1x}^{0}{}_{}{}^{2}+M_{1y}^{0}{}_{}{}^{2}}{2}}+M_{1y}^0a_+\omega _+^2M_{1x}^0a_{}\omega _{}^2`$
$`+(A_+^2+a_+^2){\displaystyle \frac{\omega _+^4}{2}}+(A_{}^2+a_{}^2){\displaystyle \frac{\omega _{}^4}{2}}].`$
Due to the averaging, the variables $`\varphi _\pm `$ are cyclic, $`P_\pm `$, the conjugate momenta of $`\varphi _\pm `$, are conserved quantities. Their expression is

$$P_\pm =\frac{L_{eff}}{\dot{\varphi }_\pm }=NA_\pm ^2\omega _\pm ^2\dot{\varphi }_\pm .$$
(54)
As the Lagrangian does not depend on the time derivatives of $`A_+`$, $`A_{}`$, $`a_+`$, and $`a_{}`$, there is no Legendre transform on these variables and the Hamiltonian reads

$`H_{eff}`$ $`=`$ $`P_+\dot{\varphi }_++P_{}\dot{\varphi }_{}+{\displaystyle \underset{i=1}{\overset{N}{}}}p_i\dot{\theta }_i^0L_{eff}`$ (56)
$`=`$ $`{\displaystyle \frac{P_+^2}{2NA_+^2\omega _+^2}}+{\displaystyle \frac{P_{}^2}{2NA_{}^2\omega _{}^2}}+{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{p_{i}^{0}{}_{}{}^{2}}{2}}`$
$`+N[{\displaystyle \frac{M_{1x}^{0}{}_{}{}^{2}+M_{1y}^{0}{}_{}{}^{2}}{2}}+M_{1y}^0a_+\omega _+^2M_{1x}^0a_{}\omega _{}^2`$
$`+(A_+^2+a_+^2){\displaystyle \frac{\omega _+^4}{2}}+(A_{}^2+a_{}^2){\displaystyle \frac{\omega _{}^4}{2}}].`$
In the absence of conjugate variables of the amplitudes $`A_\pm `$, the corresponding Hamilton’s equation are simply given, from the least action principle, by $`_{A_\pm }H_{eff}=0`$. Together with the equations (54), this leads to the following expressions for the frequencies

$`{\displaystyle \frac{d\varphi _\pm }{dt}}`$ $`=`$ $`\omega _\pm ,`$ (57)
expected of course from the previous study of the matrix $`𝐓`$, and for the amplitudes

$`A_\pm ^2`$ $`=`$ $`{\displaystyle \frac{P_\pm }{N\omega _\pm ^3}}.`$ (58)
Finally, from the equations $`_{a_\pm }H_{eff}=0`$, we find

$$a_+=\frac{M_{1y}^0}{\omega _+^2}\text{and}a_{}=\frac{M_{1x}^0}{\omega _{}^2}.$$
(59)
Reintroducing expressions (58) and (59) in the Hamiltonian (56), we end up with the *effective* Hamiltonian describing the slow motion of the particles

$`H_{eff}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{p_{i}^{0}{}_{}{}^{2}}{2}}+P_+\omega _++P_{}\omega _{}.`$ (60)
The evolution of the full system is then approximated by the dynamics of this effective Hamiltonian, the constant $`P_+`$ and $`P_{}`$ being determined by the initial conditions. However, a difficulty arises when the two eigenvalues crosses, ie when $`|M_2|`$ becomes $`0`$. Since $`\psi `$ is defined as half the phase of the complex number $`M_{2x}+iM_{2y}`$, it experiences at this point a $`\pi /2`$ jump, and consequently the eigenvectors $`X_\pm `$ are inverted: this phenomenon of eigenvalue crossing is illustrated by Fig. 9.

However, in most of the numerical experiments conducted in drh , only one fast frequency, $`\omega _{}`$ was excited. From now on, we will therefore concentrate on the case $`P_+=0`$, and will not be concerned by the phenomenon of eigenvalue crossing. Dropping the subscript for $`P`$, we consider the Hamiltonian

$`H_{eff}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{p_{i}^{0}{}_{}{}^{2}}{2}}+P\sqrt{{\displaystyle \frac{1|M_2|}{2}}}`$ (61)
The corresponding equations of motion are

$`\ddot{\theta _i}`$ $`=`$ $`{\displaystyle \frac{P}{N\sqrt{2}\sqrt{1|M_2|}}}\mathrm{sin}2\left(\theta _i+\psi \right)`$ (62)
which have to be compared with Eq. (19): this is still a pendulum equation, but the potential amplitude, depending now self consistently on the particles motion through $`M_2`$, explains the presence of untrapped particles, shown in Fig. (5b). In addition, we have now the right expression for the frequency of the fast oscillations $`\omega =\omega _{}`$, and a new conserved quantity has been identified: the adiabatic invariant $`P`$. Using Eq. (58), we obtain the $`A_{}(M_2)`$ relation which, together with $`\omega (M_2)`$, is perfectly verified by numerical simulation as attested by Fig. 10.

The full efficiency of this procedure is that it preserves completely the Hamiltonian structure of the problem, making it well suited for a statistical treatment, presented below. However, it is important to emphasize before, that all fast variables having disappeared, this procedure induces a huge gain (of order $`1/\epsilon `$) for numerical simulations. This allows to study accurately the statistics and the dynamics of this effective Hamiltonian up to extremely long times, giving conclusion directly relevant to the original model.

### 4.3 Study of the effective Hamiltonian

#### 4.3.1 Statistical mechanics of the effective Hamiltonian

We would like to describe the bicluster as an equilibrium state of the effective Hamiltonian (61). Fortunately, due to the special form of the potential, which depends only on the global variable $`|M_2|`$, the statistical mechanics of the effective Hamiltonian are exactly tractable in the microcanonical ensemble.

There are two conserved quantities: the total energy $`E`$ and the total angular momentum. As the latter only creates a global rotation of the system totally decoupled from the rest of the dynamics, we can consider the total energy without restriction. The density of states $`\mathrm{\Omega }(E)`$ at a given energy, written as an integral in which the energy is split in two parts, a kinetic and a potential one, has the following expression

$`\mathrm{\Omega }(E)`$ $`=`$ $`{\displaystyle 𝑑V\mathrm{\Omega }_{kin}(EV)\mathrm{\Omega }_{pot}(V)}`$ (63)
$``$ $`{\displaystyle d|M_2|\mathrm{\Omega }_{kin}\left(EV(|M_2|)\right)\mathrm{\Omega }_{conf}(|M_2|)}`$
where the sign $``$ accounts for an unimportant constant. $`\mathrm{\Omega }_{kin}`$ (resp. $`\mathrm{\Omega }_{pot}`$) is the density of states corresponding to the kinetic (resp. potential) part of the Hamiltonian, and $`\mathrm{\Omega }_{conf}`$ is the density of angular configuration giving rise to a given $`|M_2|`$: since the potential only depends on $`|M_2|`$, $`\mathrm{\Omega }_{pot}(V)`$ is directly proportional to $`\mathrm{\Omega }_{conf}(|M_2|)`$. Their expressions are

$`\mathrm{\Omega }_{kin}(E)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \underset{i}{}}dp_i\delta \left({\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}p_i^2E\right)`$ (64)
$`\mathrm{\Omega }_{conf}(|M_2|)`$ $`=`$ $`{\displaystyle _0^{2\pi }}{\displaystyle \underset{i}{}}d\theta _i\delta ({\displaystyle \frac{1}{N}}\left({\displaystyle \underset{i}{}}\mathrm{cos}2\theta _i\right)^2`$ (65)
$`+{\displaystyle \frac{1}{N}}\left({\displaystyle \underset{i}{}}\mathrm{sin}2\theta _i\right)^2N|M_2|^2).`$
Using the classical result for a perfect gas, we obtain

$`\mathrm{\Omega }_{kin}(E)`$ $`=`$ $`{\displaystyle \frac{2\pi ^{N/2}}{\mathrm{\Gamma }(N/2+1)}}E^{N/2}.`$ (66)
To compute $`\mathrm{\Omega }_{conf}(|M_2|)`$, it is possible to use the inverse Laplace transform of the Dirac delta function

$$\delta (x)=_\mathrm{\Gamma }𝑑pe^{px},$$

where $`\mathrm{\Gamma }`$ is a path in the complex plane running from $`\mathrm{}(p)=\mathrm{}`$ to $`\mathrm{}(p)=+\mathrm{}`$. The expression of $`\mathrm{\Omega }_{conf}`$ becomes

$`\mathrm{\Omega }_{conf}(|M_2|)`$ $`=`$ $`{\displaystyle _\mathrm{\Gamma }}𝑑pe^{Np|M_2|^2}{\displaystyle _0^{2\pi }}{\displaystyle \underset{i}{}}d\theta _ie^{NpM_2^2}.`$ (67)
To compute the integral over the angles, we use now the Gaussian (also called Hubbard-Stratanovich) transform

$`e^{Npm^2}`$ $``$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑ue^{N\frac{u^2}{4p}+Nmu}.`$
This decouples the integration over the angles $`\theta _i`$, and we have

$`\mathrm{\Omega }_{conf}(|M_2|)`$ $``$ $`{\displaystyle _\mathrm{\Gamma }}𝑑pe^{Np|M_2|^2}{\displaystyle 𝑑u𝑑ve^{\frac{N}{4p}(u^2+v^2)}}`$ (68)
$`{\displaystyle _0^{2\pi }}{\displaystyle \underset{i}{}}d\theta _ie^{N_iu\mathrm{cos}2\theta _i+v\mathrm{sin}2\theta _i}`$
$``$ $`{\displaystyle _\mathrm{\Gamma }}𝑑pe^{Np|M_2|^2}{\displaystyle _0^{\mathrm{}}}r𝑑re^{\frac{Nr^2}{4p}}I_0(r)^N`$ (69)
where $`r`$ is the polar radius associated with $`u`$ and $`v`$, and $`I_0`$ is the modified Bessel function of order $`0`$. Using the saddle point method on $`p`$ and $`r`$, we have finally

$`\mathrm{\Omega }_{conf}(|M_2|)`$ $``$ $`e^{N\left(p^{}|M_2|^2+\frac{r_{}^{}{}_{}{}^{2}}{4p^{}}\mathrm{ln}I_0(r^{})\right)}`$ (70)
$`=`$ $`e^{N\left(2p^{}|M_2|^2\mathrm{ln}I_0(r^{})\right)}`$ (71)
where $`p^{}`$ and $`r^{}`$ are implicitly defined by

$`{\displaystyle \frac{r^{}}{2p^{}}}`$ $`=`$ $`{\displaystyle \frac{d\mathrm{ln}I_0(p)}{dp}}_{|r^{}}`$ (72)
$`{\displaystyle \frac{r_{}^{}{}_{}{}^{2}}{4p_{}^{}{}_{}{}^{2}}}`$ $`=`$ $`|M_2|^2.`$ (73)
Using this result, Eq. (63) may be evaluated once again by the saddle point method, which leads after some simple algebra to the equilibrium value $`|M_2|^{}`$, defined by the following equation

$`r^{}(|M_2|^{})`$ $`=`$ $`{\displaystyle \frac{1}{4\sqrt{1|M_2|^{}}\left({\displaystyle \frac{\sqrt{2}E}{P}}\sqrt{1|M_2|^{}}\right)}}.`$ (74)
As shown by Fig. 12, this equation has only one solution for any ratio $`E/P`$; this indicates an absence of phase transition. As this solution is non zero, a biclusterization will always take place; nevertheless in the large energy limit the particles are almost free and $`|M_2|^{}`$ goes of course to zero.

Once $`|M_2|^{}`$ is known, it is easy to complete the description of the equilibrium state. The temperature $`T`$ is given by

$`T`$ $`=`$ $`{\displaystyle \frac{2}{N}}E_c={\displaystyle \frac{2}{N}}EP\sqrt{{\displaystyle \frac{1|M_2|^{}}{2}}}`$ (75)
the distribution of velocities is a Maxwellian with temperature $`T`$, and the distribution of angles has a gibbsian shape $`\rho (\theta )e^{V(\theta )/T}`$ with the potential

$`V(\theta )`$ $`=`$ $`{\displaystyle \frac{P}{N2\sqrt{2}\sqrt{1|M_2|^{}}}}\left(1\mathrm{cos}(2\theta +2\psi )\right)`$ (76)
This potential is inferred from the equations of motion (62).

#### 4.3.2 Relaxation to equilibrium

We have therefore now a complete description of the statistical equilibrium states of the effective Hamiltonian, governed by long range interactions. It has been noticed by various authors that the dynamics of such systems may sustain long lived metastable states before relaxing to equilibrium Antoni ; AnteneodoTsallis ; Latora1 ; gouda . Before comparing the ”effective equilibrium” with the structure created by the dynamics of the real Hamiltonian, it is thus necessary to study the relaxation to equilibrium. As the relaxation properties of long range interacting systems is in itself an important problem, we will consider them now. We will in addition illustrate which statistical properties of the equilibrium distribution is expected to be observed in the real dynamics, for large $`N`$ and on time scale reasonable for numerical computations.

Fig. 13 illustrates the approach to equilibrium of the effective Hamiltonian, for three different particle numbers, with initially immobile and homogeneously distributed particles (this corresponds to the typical initial condition used in the real Hamiltonian in Antoni ; drh  for instance). We have represented $`<M_{2n}>(\tau )`$, the time averages of $`M_{2n}`$ from initial time $`0`$ to time $`\tau `$, for $`n`$ equal to 2, 4 and 6. Temporal fluctuations are thus not visible. We first observe finite size effects concerning the equilibrium values : for instance $`M_6`$ converges for the small system $`N=200`$ to a value larger than $`M_6^{}`$, the equilibrium value.

More interestingly, we observe also that whereas $`|M_2|`$ quickly reaches its equilibrium value, the relaxation time of the successive moments $`|M_4|`$, $`|M_6|`$ strongly depends on the system size, and presumably diverges with $`N`$. The actual dependence of the relaxation time of each moment, on the particle number, will not be studied in this paper. It would for instance address the issue of whether non equilibrium distributions may be observed in the thermodynamic limit. Such phenomena have indeed already been observed in other long range interacting systems Latora1 ; Latora2 .

We conclude that for a large, but finite, particle number, the moments of the distribution converge toward their equilibrium values. When $`N`$ increases, the relaxation time for the high order moments may however be larger than computationally achievable times and, we thus expect the simulation to exhibit very slowly evolving non equilibrium structures of the effective Hamiltonian. Authors of drh  have indeed found a distribution of particles different from the Gibbsian shaped distribution (76) predicted by the equilibrium thermodynamics of the effective Hamiltonian (see Fig. 8).

On the contrary our analysis shows that $`|M_2|`$ quickly reaches its equilibrium value $`|M_2|^{}`$. For this reason, given the computational time achievable, in the next section, we only use this statistical equilibrium indicator for the study of the structure. This will lead to the correct prediction of the equilibrium repartition of potential and kinetic energy, and of its dependance from initial condition.

### 4.4 Comparison with the full Hamiltonian

If both time scales are clearly separated, Fig. 11 emphasizes the striking agreement between the effective and the real dynamics, for short times. We will show in this section that the effective Hamiltonian also provides a good description for the long time dynamics.

For this purpose, we will first compare the evolution of the non-equilibrium angular density distribution $`\rho (\theta )`$ for both dynamics. The results are reported on Fig. 14 were circles show the repartition given by the original Hamiltonian (1) whereas the dashed line shows the results of the effective dynamics (61). This picture clearly emphasizes that the non equilibrium character of the dynamics is fully described by the effective dynamics, and we can easily conclude that the effective Hamiltonian reflects very well the dynamics of the real Hamiltonian.

As discussed in the previous section, according to the computational times used, only the second momentum of the distribution $`M_2`$ should correspond to the equilibrium one. Let us therefore analyze the dependence of the equilibrium value of $`M_2`$ on the initial conditions. The typically used initial condition for the numerical simulations in the original system (1) is an initially quasi homogeneous distribution of particles with zero velocity drh . It corresponds to the ratio $`e/\sqrt{2}P=1`$ for the effective Hamiltonian. According to Fig. 12 , this implies an equilibrium value $`|M_2|^{}0.51`$ in perfect agreement with the numerical result reported earlier drh .

From these results, it is also possible to explain the caloric curve $`T1.3e`$, reported earlier drh . The total energy of the full system is divided in three parts: the potential energy of the small oscillations, the kinetic energy of the small oscillations and the kinetic energy of the slow movements. The two first parts are equal in average and form the potential part of the effective Hamiltonian, whereas the last one is the kinetic energy of the effective Hamiltonian. Using this remark together with the values of $`e/P`$ and $`|M_2|^{}`$ at equilibrium, it is easy to derive the temperature/energy relationship, using (75). In the case investigated in drh , we find precisely $`T=1.3e`$.

We thus conclude that the dynamics of the effective Hamiltonian parametrizes very well the dynamics of the real Hamiltonian, for short as well as for long time. This allows us to predict statistical properties of the initial system, as for instance the asymptotic value of $`M_2`$ or the partitioning between kinetic and potential energy. Moreover, let us recall that the effective Hamiltonian gives the opportunity to study numerically the relaxation towards equilibrium of the bicluster, whereas it was not possible in the real dynamics, because of computational limitations (let us note that the ratio of the typical time scale of the two dynamics is of order 100 or larger).

All these findings are great successes of this approach. Let us nevertheless comment some points that we have not yet addressed.

The first one concerns the limit of validity of our multiple time scale analysis. It should be noticed that not all initial conditions with small energy would lead to the formation of the bicluster. From our analysis, we can conclude that the class of initial conditions leading to this formation is the one compatible with the ansatz used: however a precise description of this class is not known. In particular, we have not study the threshold value of $`\epsilon `$ beyond which such a description breaks down. For $`\epsilon `$ going toward one, the ansatz (41) will first lead to a nonlinear evolution of the slow variable $`f_i`$ in place of the linear one given by Eq. (46). We think that the critical value $`\epsilon _c`$ above which this nonlinear dynamics do not have any more periodic solution should correspond to the critical value above which the bicluster can not form, explaining the transition observed in drh . Such an hypothesis should be tested.

A second point is the phenomenon of level crossing discussed in section 4.2, leading to an exchange of excitation of the two modes of the system. This phenomenon is very similar to level crossing in quantum mechanics adiabatic problems. It represents a resonance, local in time, in which an interaction between the various modes may occur, leading to a modification of the adiabatic invariant. For sake of simplicity, we have concentrated our work on the case of excitation of only the smallest frequency, so that such level crossing could not occur. A more general study would however be of interest.

The last point concerns the validity of our approach for infinite times. As noted in the introduction, a recent paper FFR  has shown a long-time degradation of the bicluster, for a very small number of particles, suggesting its transient non-equilibrium nature in that case. For a reasonable number of particles, our numerical results show that such a degradation does not exist, for computationally achievable times. This point is thus not addressable numerically. From a theoretical point of view, results on much more simple systems, show that adiabatic invariants as described here, should be conserved for a very long time (exponential when $`\epsilon `$ goes to 0). Moreover, such very long time stability results, would not give any hint on what happens for larger times (stability or instability).

## 5 Conclusion

The surprising formation and stabilization drh  of the bicluster in the HMF dynamics, in contradiction with statistical mechanics prediction, is now understood: the small collective oscillations of the bulk of particles create an effective double-well potential, in which the particles evolve.

On a first stage the dynamics can be described in the context of a forced Burgers’ equation. The particles trajectories caustics form infinite density regions explaining an initially diverging density. Early times dynamics is therefore very similar to the structure formation in Eulerian coordinates for a one dimensional self-gravitating system noullez . As in this case, because of the confining potential, the particles can not move apart as easily after the first caustic has been formed, as in the case of the free motion on the plane noullez .

On much larger timescales, in order to parametrize the fast oscillations of the bulk particles, we have performed a variational multiscale analysis. We then have obtained an Hamiltonian effective model describing the particle slow motion. This description is in very good agreement with the initial dynamics. The effective dynamics allows numerical simulations on time scale much larger than the initial dynamics, as the rapid oscillations have been filtered out. We have performed the statistical mechanics of this new system. The results thus give a statistical explanation of the bicluster formation and stabilization. The equilibrium properties can then explain the mean value of the second moment of the distribution and the repartition between potential and kinetic energy. Even if this analysis has only been sketched in this work, the effective dynamics is also a powerful tool to study the very slow relaxation process versus the real equilibrium structure.

Many physical systems share the main properties of this HMF dynamics : very fast oscillations self interacting with a slower motion. In addition to the plasma problem already cited and always in the context of long range interacting systems, we may for instance cite the problem of interaction of fast inertia gravity waves with the vortical motion, for the rotating Shallow Water or the primitive equation dynamics, in the limit of a small Rossby Number embidmajda . The main interest of our study is to provide a toy model in which such a complex dynamics can be treated and analyzed extensively using powerful theoretical tools. We point out in particular the usefulness of a variational approach in the procedure of averaging the rapid oscillations. This toy model also permits a clear view of the usual problems of such dynamics
(i) Resonances (or level crossing) in averaging procedures
(ii) Unusual relaxation processes due to the long range nature of the interactions.
The HMF model is to our knowledge the simplest one in which such phenomena occur. The study of these phenomena is thus a natural extension of our work and should be of interest in analyzing more complex systems.

###### Acknowledgements.

We would like to thank M-C. Firpo, P. Holdsworth, S. Lepri, J. Lebowitz, F. Leyvraz, A. Noullez, for very helpful advices. This work has been partially supported by EU contract No. HPRN-CT-1999-00163 (LOCNET network), the French Ministère de la Recherche grant ACI jeune chercheur-2001 N 21-31, and the Région Rhône-Alpes for the fellowship N 01-009261-01. This work is also part of the contract COFIN00 on Chaos and localization in classical and quantum mechanics.

## Appendix A: Variationnal description of a simple adiabatic dynamics

Let us consider a slowly modulated harmonic oscillator with Lagrangian

$$L(\dot{\theta },\theta ,t,\tau )=\frac{\dot{\theta }^2}{2}\omega ^2(\tau )\frac{\theta ^2}{2}$$
(77)
where $`\tau =\epsilon t`$. We consider the following ansatz

$$\theta (t,\tau )=A(\tau )\mathrm{sin}\phi (t)$$
(78)
with $`\dot{\phi }=O(1)`$ and $`\dot{A}=\epsilon dA/d\tau `$, but, contrary to the usual asymptotic expansion on the equation of motion, we will consider the variational approach proposed by Witham witham . The idea is to separate the two different time scales of the motion at the level of the action. The average of the fast variable $`t`$ yields the following effective Lagrangian:

$$=L_t=\frac{A^2}{4}\left(\dot{\phi }^2\omega ^2\right).$$
(79)
It is thus straightforward to obtain the equation of motion of the effective Lagrangian $`(A,\dot{A},\phi ,\dot{\phi })`$. They read

$`{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{}{\dot{A}}}\right)`$ $`=`$ $`{\displaystyle \frac{}{A}}`$ (80)
$`{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{}{\dot{\phi }}}\right)`$ $`=`$ $`{\displaystyle \frac{}{\phi }}`$ (81)
i.e.

$`0=A\left(\dot{\phi }^2\omega ^2\right)`$ $``$ $`\dot{\phi }=\omega `$ (82)
$`{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{A^2\dot{\phi }}{2}}\right)=0`$ $``$ $`\dot{A}^2\dot{\phi }=\text{cste}`$ (83)
This method emphasizes a new constant of the motion $`A^2\dot{\phi }`$, called adiabatic invariant, that is not exhibited by the usual asymptotic expansion on the equations of motion corresponding to the original Lagrangian (79). This may be of great interest in complex systems. Moreover, this invariant necessarily appears as the angle $`\varphi `$ associated to the fast variable $`t`$ is always cyclic after averaging. The same reasoning also apply for larger order ansatzs, showing that an adiabatic invariant should exist at any order in $`\epsilon `$.

In addition, this method preserves the Hamiltonian character of the problem allowing for instance a statistical mechanics treatment as presented in this paper.

## Appendix B: Eigenvalues and eigenvectors of the matrix $`𝐓`$

To solve the linear second order equation (46), we need to diagonalize the $`N`$th-order circulant matrix $`𝐓`$ defined by $`T_{ij}=\mathrm{cos}(\theta _i\theta _j)/N`$. We introduce the two vectors $`𝐗(\psi )`$ and $`𝐘(\psi )`$, with coordinates $`X_i(\psi )=\mathrm{cos}(\theta _i+\psi )`$ and $`Y_i(\psi )=\mathrm{sin}(\theta _i+\psi )`$, and $`\psi `$ an arbitrary phase. $`𝐓`$ is therefore the sum of the projections along these vectors

$$𝐓=\frac{1}{N}\left(𝐗(\psi )^t𝐗(\psi )+𝐘(\psi )^t𝐘(\psi )\right).$$
(84)
This proves that the image of $`𝐓`$ is of dimension two, and that only two eigenvalues $`\lambda _\pm `$ are nonzero. Let us restrict now to the two dimensional problem in the $`Vect(𝐗,𝐘)`$ plane (this plane does not depend on $`\psi `$ of course), and let us choose $`\psi `$ such that $`X`$ and $`Y`$ are orthogonal, i.e.

$$(𝐗.𝐘)=\underset{i}{}\mathrm{cos}(\theta _i+\psi )\mathrm{sin}(\theta _i+\psi )=0.$$
(85)
This definition of $`\psi `$ leads to

$$|M_2|=\frac{1}{N}\underset{i}{}\mathrm{cos}(2\theta _i+2\psi ).$$
(86)
Once we have the two eigenvectors $`𝐗_+=𝐗(\psi )`$ and $`𝐗_{}=𝐘(\psi )`$, their associated eigenvalues are defined by the following relationship, valid $`i`$:

$$\underset{j}{}\frac{1}{N}\mathrm{cos}(\theta _i\theta _j)X_{\pm i}=\lambda _\pm X_{\pm i}.$$
(87)
leading directly to

$$\lambda _\pm =\frac{1\pm |M_2|}{2}.$$
(88)


format: text

text:
# Hubble Space Telescope Imaging of the Post-Starburst Quasar UN J1025-0040: Evidence for Recent Star Formation

## 1 Introduction

Exciting progress has been made bolstering our understanding of the connection between active galactic nuclei (AGNs) and starbursts. It appears that essentially every large galaxy harbors a supermassive dark object (presumably a black hole) some 0.15% of the mass of its spheroidal component (Magorrian et al. 1998; Gebhardt et al. 2000a), although the black hole mass perhaps correlates more tightly with the velocity dispersion of the host bulge (e.g., Gebhardt et al. 2000b; Ferrarese & Merritt 2000a). These are plausibly “dead” quasars, which once burned while undergoing high accretion episodes. Furthermore, black hole masses for “live” quasars determined from reverberation techniques also correlate with the bulge mass and velocity dispersion (Gebhardt et al. 2000b; Ferrarese et al. 2001). Consistent with these results, there exists an upper limit to the luminosity of AGNs that depends on their host galaxy stellar mass (e.g., McLeod, Rieke, & Storrie-Lombardi 1999). Finally, there is the global similarity in the strong evolution of the quasar luminosity function and the star-formation rate (e.g., Boyle & Terlevich 1998).

These points have sparked much recent theoretical interest in hierarchical galaxy formation and related scenarios featuring direct relationships between star formation and the growth/fueling of central black holes (e.g., Burket & Silk 2001, and references therein). AGN activity, in fact, may reflect a fundamental stage of galaxy evolution.

The discovery of a large population of ultraluminous infrared galaxies (ULIRGs; $`L_{IR}`$ $`>10^{12}L_{}`$) that are strongly interacting merger systems heated by both starburst and AGN power sources led Sanders et al. (1988) to hypothesize that ULIRGs represent the initial dust-enshrouded stage of a quasar (see also Sanders & Mirabel 1996 for a review). The evidence for this evolutionary scenario is circumstantial, but the hypothesis is supported by HST images of warm ULIRGs and optical QSOs that show similarities such as knots of recent star formation (e.g., Surace et al. 1998; Boyce et al. 1996; Bahcall et al. 1997). But Boyce et al. (1996) suggested that the $`IRAS`$-selected quasars are the products of young spiral mergers, not old mergers. Recent $`ISO`$ results may also be incompatible with a strict evolutionary scenario: Genzel et al. (1998) surveyed ULIRGs and found that half contain simultaneously an AGN and starburst activity in a 1-2 kpc circumnuclear ring (these very young star-forming regions have ages 10-100 Myrs), but that the presence of an AGN is not correlated with the merger age (see also Canalizo & Stockton 2001).

During the course of identifying quasar candidates in the AAT-2dF quasar survey (Croom et al. 2001) coincident with radio sources in the NRAO VLA Sky Survey (NVSS; Condon et al. 1998), we discovered UN J1025$``$0040 (Brotherton et al. 1999a). This object’s unusual spectrum shows the broad Mg II emission line of a quasar as well as the large Balmer jump and deep Balmer absorption lines reminiscent of an A-type star, all at $`z=0.634`$. The spectrum may be categorized as a “Q + A” (quasar plus post-starburst) spectrum similar to those of post-starburst galaxies which have “k + A” spectra due to the combination of young and old stellar populations. We modeled the stellar content of UN J1025$``$0040 using stellar synthesis population models (Bruzual & Charlot 1996), and found that a 400-Myr-old instantaneous starburst provides an excellent match to the Balmer jump and stellar absorption lines. In conjunction with the broad emission lines and power-law spectrum of a blue quasar, the observed spectrum could be completely accounted for, although not uniquely. A population of very young stars could also be present, but such a population lacks strong features in the observed spectral region to betray its presence (ultraviolet spectroscopy or imaging analysis, as in this paper, is needed).

The estimate of the initial mass of the 400-Myr-old starburst was nearly $`10^{11}`$ $`M_{\mathrm{}}`$, or the mass of our own Milky Way Galaxy. Subtracting off this model left a residual spectrum that resembles a normal if rather excessively blue quasar. The total energy budget of both components appeared similar, $`10^{11.6}`$ $`L_{\mathrm{}}`$. Running the stellar synthesis population model in reverse, when the starburst was only 50 Myr old it would have had a bolometric luminosity of $`10^{12.4}`$ $`L_{\mathrm{}}`$ and qualified as an ULIRG.

Deep Keck Ks-band imaging with 0.5″ seeing failed to resolve quasar and stellar components although it did reveal a companion galaxy 4.2″ away that turned out to have a dominant 800-Myr-old stellar population (Canalizo et al. 2000). Faint \[O II\] $`\lambda `$3727 knots are also present between UN J1025$``$0040 and the companion. Certainly this is an interacting system, as are essentially all ULIRGs.

With these results in hand, we were awarded HST cycle 9 time to image UN J1025$``$0040 with WFPC2 and attempt to resolve the starburst from the quasar. The sub-pixel dithered resolution of the WFPC2 PC achieves a spatial resolution of about 600 parsecs FWHM at $`z=0.634`$ ($`H_0=75`$ km s-1Mpc-1, $`q_0=0`$), and we anticipated seeing a $``$kpc-scale starburst ring or a double nucleus.

We do not resolve any such structures (§ 2). We marginally resolve extended emission (stars) around a point source (the AGN) and see a disturbed morphology suggestive of a nuclear merger remnant (§ 3). The ratio of the point-to-extended light differs from the spectral analysis of Brotherton et al. (1999a) and indicates the presence of a previously unsuspected young or on-going starburst (§ 4). We discuss the implications of these results (§ 6) and summarize our conclusions (§ 7).

## 2 HST Images

We imaged UN J1025$``$0040 in two broad bands, F814W in the red longward of the redshifted Balmer jump, and F450W in the blue shortward of the redshifted Balmer jump. UN J1025$``$0040 was centered on the Planetary Camera (PC) detector of WFPC2 (0.046 arcsec per pixel), and multiple images at each band were taken using a sub-pixel dithered box pattern for maximum resolution. Eight 500 sec observations were made on 13 Apr 2001 (UT) through the F814W filter, and sixteen 500 sec observations on 20 Apr 2001 (UT) through the F450W filter.

We began with the pipeline-calibrated images which already had undergone bias subtraction and flat fielding. We combined our dithered images using two different methods, the DRIZZLE algorithm of Fruchter et al. (1997), and the algorithm of Lauer (1999). Both algorithms produce images with twice the original number of pixels along each axis, 0.023 arcsec per pixel in the case of the PC. These combination procedures also efficiently clean cosmic rays. An examination of the profiles of the resulting images showed that the algorithm of Lauer (1999) produced sharper images. For the blue F450W images, the FWHMs of cross-cuts through the point-source-dominated peaks were 3.165 pixels versus 3.875 pixels, and for the red F814W image the FWHMs were 4.537 and 5.042 pixels, respectively, with the Lauer reduction superior in each case. We therefore used the images produced by the Lauer (1999) algorithm for our analyses, which are shown in Figures 1 and 2.

The compact companion galaxy 4.2″ to the south-southwest is detected and resolved for the the first time in these HST images. In the F450W image the companion spans a total spatial extent of 0.23$`\mathrm{}`$ and in the F814W image 0.35$`\mathrm{}`$; the size in the red image is then 2 kpc for $`h=0.75`$, which was approximately the limit reported by Canalizo et al. (2000). Table 1 gives the photometry of the integrated flux of UN J1025$``$0040 and its companion through the F450W and F814W filters. Figure 3 shows that the photometry is in good agreement with that reported by Canalizo et al. (2000), indicating that the AGN has not varied significantly between the two epochs. This point will be important in reaching some of our conclusions below. The magnitudes of the companion galaxy are also entirely consistent with our previous ground-based photometry.

## 3 Image Analysis

From the fully reduced images it is apparent that we have failed to clearly resolve any significant starburst structures (e.g., a second nucleus or a ring). The images do resolve much of the centrally concentrated starburst, which is extended over a diameter of $``$1 kiloparsec around the sharper point source nucleus that is the AGN. Some spurs of low surface brightness extended outward, the most notable being toward the northeast. A low surface brightness plateau extends to the west about 0.5 arcsec. The outermost surface brightness contour in the red image shows is perhaps suggestive of a merger.

There were no suitable stars in the image field of the PC, so we used the Tiny Tim software (Krist & Hook 1999) to produce model PSFs in conjunction with the spectral shape of the point source, a blue quasar, in this case using a composite quasar spectrum (Brotherton et al. 2001). We tried different values of the jitter parameter to find the best estimate for the PSF in each image, as telescope breathing and detector effects make the PSF impossible to model perfectly. The jitter added was 14 mas for the F450W model, and 18 mas for the F814W model.

Once we had good estimates of the PSFs, we used the Richard-Lucy algorithm to deconvolve our images. These deconvolved images did not reveal any significant additional detail. We have a disturbed nuclear region with a centrally concentrated, amorphous structure marginally resolved at the PC scale. The very low level extended (arcsecond scale) starlight is energetically negligible compared to the centrally concentrated light. This surrounding “fuzz” has colors that are very red, $`m(F450W)m(F814W)=2.0`$ mag in the northeast and $`m(F450W)m(F814W)=1.6`$ mag in the west, indicating an old (1+ Gyr) stellar population.

We also created PSF-subtracted images to determine the ratio of point source to extended flux. The minimum amount of point source present was estimated by subtracting off increasing amounts of the PSF model until a smooth residual was obtained that did not have a central dip. The maximum amount was estimated by subtracting off increasing amounts of the PSF model until zero flux was reached in the residual. This uncertainty dominates those present in the choice of PSF model. We estimate that the extended emission represents $`52\pm 4\%`$ of the total flux in the F450W image, and $`69_5^{+4}\%`$ in the F814W image.

## 4 Spectral Analyses: Quasar and Starbursts

Figure 4 shows an adaptation of a figure from Brotherton et al. (1999a), in which a 400-Myr-old instantaneous starburst accounts for the stellar light in UN J1025$``$0040 and a rather blue AGN accounts for the remainder. This previous analysis indicated that the starburst contributes some two-thirds of the red light, and one quarter to one third of the blue light. The size of the Balmer jump and the depths and equivalent width ratios of the stellar absorption lines provided excellent constraints in the red part of the spectrum, but there are few stellar features in the blue to check the model.

The image analysis described in the previous section provides an independent estimate of the fractions of AGN and starburst in each band. While there is good agreement in the red F814W band between the starburst fractions determined by the image analysis and the previous spectral analysis, the blue F450W band shows a surprise. The starburst fraction in the blue is half the total, not $``$30% as expected, a significant difference.

We have an excess of extended blue light, so it cannot arise from direct AGN light. The polarization is low (1-2%; Brotherton et al. 1999a) so it cannot arise from indirect, scattered light from the AGN, and the size of the scattering region would have to be implausibly large.

Single-age starbursts cannot simultaneously explain the spectral features in the red part of the spectrum and the blue-red spectral energy distribution, even when allowing the several free parameters to vary (e.g., metallicity, dust, initial mass function).

Very young stars can contribute significant blue light, while at the same time having only weak Balmer jumps and weak absorption lines. Such a population, added to a 400-Myr-old population, can make a composite starburst spectrum bluer yet allow it to retain essentially the same Balmer jump and absorption line equivalent width ratios (which are unaffected by the diluting quasar light). Figure 5 illustrates again the spectrum of UN J1025$``$0040 as in Figure 4, but also shows the total light photometry, the estimated starburst photometry, and a composite starburst population model. The composite model was created using Bruzual & Charlot (1996) stellar synthesis population models with two instantaneous starbursts (Salpeter IMF, solar metallicity, no dust), one 400 Myr old, the other 40 Myr old or younger (10 Myr old in the case of Fig. 5), scaled to match the combination of imaging and spectral constraints. Young starbursts older than 40 Myr are increasingly red, have larger Balmer jumps, and more pronounced absorption lines; such populations are still bluer than the 400 Myr old model and capable of meeting the imaging constraints when added in the correct proportions, but will fail to simultaneously match the spectral constraints of the Balmer jump and absorption lines.

As discussed previously by Brotherton et al. (1999a), varying the model parameters (primarily metallicity) can result in somewhat different ages. In the face of increasing complexity of this system and the limited data so far available, determining the exact parameters and ages is not yet possible. Nevertheless, we have a robust result that a significantly younger stellar population is present in the nucleus of UN J1025$``$0040.

The residual quasar spectrum implied by this analysis is not as blue as originally suggested by Brotherton et al. (1999a). The quasar continuum can be approximated by a power-law continuum where the optical spectral index $`\alpha =0.5`$, where $`f_\nu \nu ^\alpha `$. This is a very typical optical power-law index (e.g., Francis 1993).

How young (and hence also how massive) is the young starburst? How does its presence affect our previous conclusions regarding the older, more massive starburst? Figure 6 shows a range of viable models. All models shown are combinations of dust-free instantaneous starbursts employing Salpeter initial mass functions and solar metallicities, and the composite quasar spectrum of Brotherton et al. (2001) based on the FIRST Bright Quasar Survey (FBQS; White et al. 2001). The combination of bright optical selection and deep radio selection in the FBQS is rather similar to that used for the 2dF-NVSS selection that identified UN J1025$``$0040. This composite quasar spectrum has a power-law index of $`\alpha =0.46`$, consistent with that determined above.

The scaling, and hence mass, of the 400-Myr-old starburst is rather insensitive to the age of the younger population. For the updated photometry here and for $`H_0=75`$ km s-1Mpc-1, $`q_0=0`$, the oldest composite starburst model of Figure 6 contains a 400-Myr-old starburst with $`5.7\times 10^{10}M_{}`$, plus a 40-Myr-old starburst with $`6.7\times 10^9M_{}`$, or 10% of the mass of the older starburst. For the youngest models of Figure 6 (0 and 1.5-Myr-old instantaneous bursts), the starburst masses are down another order of magnitude, or about 1% of the mass of the older starburst. At a good fraction of 10$`{}_{}{}^{9}M_{}`$, this still qualifies as a very massive starburst in its own right. The exact parameters are not well determined and require additional constraints in the form of ultraviolet spectroscopy. The caveats of Brotherton et al. (1999a) regarding the uncertainties on the age and mass of the older starburst should also be kept in mind.

## 5 Discussion

The primary aim of obtaining the HST images was to resolve the starburst from the AGN and identify any morphological structures that might provide a guide as to the origin of such an extremely massive starburst. In that, we have been partially successful, insofar as we can distinguish a point source from a more diffuse region of centrally concentrated starlight. Since the ratio of starlight to quasar in the red part of spectrum is just as predicted from the spectral analysis of Brotherton et al. (1999a), we are probably not missing any core of starlight still unresolved from the AGN.

We were not successful in identifying any clearly recognizable star-forming structures, such as a starburst ring as has been seen in some other AGNs (e.g., NGC 985, Rodriguez Espinosa & Stanga 1990). At these redshifts the HST resolution is only probing the marginally sub-kiloparsec scale, so it is still possible that such an organized structure exists on smaller scales but is not yet resolvable as such.

These shallow images do not reveal as much detail outside the nucleus as do the deeper low-resolution ground-based images (Brotherton et al. 1999a; Canalizo et al. 2000). The nuclear structure is revealed to be non-uniform and suggestive of an evolved merger product, perhaps similar to Arp 220 (e.g., Scoville et al. 1998), but not as dusty. Whether or not UN J1025$``$0040 is an actual merger product or merely disturbed and interacting will likely require much deeper high-spatial-resolution images. Interpreting UN J1025$``$0040 as a merger product, however, would be entirely consistent with the evolutionary hypothesis of Sanders et al. (1988) in which the end result of massive mergers would be such a post-starburst quasar. This evolutionary hypothesis has been recently boosted by detailed work on the IRAS 1-Jy sample of ULIRGs (e.g., Veilleux, Kim, & Sanders 1999) which shows an increasing fraction of luminous AGNs with increasing merger age. This is a large, complete sample that includes primarily the most luminous ULIRGs for which the evolutionary hypothesis may hold true overall; contrary results, like those of Genzel et al. (1998), may result from small samples including less luminous ULIRGs. This is an ongoing area of debate.

The new result of the presence of a young ($``$40 Myr) starburst adds to this debate, indicating that specific systems can indeed be complex. Canalizo & Stockton (1997) suggested that in an evolutionary progression, the age of a starburst is a clock indicating the age of the AGN activity (Canalizo & Stockton 2001 strengthens the original arguments but also provides additional caveats). This assumes that the fueling of the nuclear starburst and the central black hole occurs at the same time, and that the black hole continues to be fueled. Determining the ages of starbursts in AGN systems would then permit the quasar lifetime to be estimated, and the quasar duty cycle as well. These numbers would complement those determined from combining the evolution of the quasar luminosity function and the population of relic black holes today in the local universe. The complication of our new result is that there has been a substantial starburst in UN J1025$``$0040 much more recently than 400 Myr, and this starburst may represent the appropriate clock for the current activity. The AGN may well not shine steadily, but undergo bursts that correspond to events that provide fuel and also incite star formation.

Our measurements of the starburst masses in the nucleus of UN J1025$``$0040, a total of 6-7 $`\times `$ 1010 M, represent an estimate or at least lower limit to the host galaxy bulge mass. We can estimate the AGN properties, black hole mass and fraction of the Eddington luminosity at which it shines, using techniques alluded to in the introduction (following the formalism and choice of parameters of Merritt & Ferrarese 2001, which includes our chosen cosmology). A black hole mass is estimated by assuming Keplerian motions of ionized gas in a spherical broad emission line region (BLR), with a velocity $`v_{BLR}=\sqrt{3}/2\mathrm{FWHM}(H\beta )`$ at a distance $`R_{BLR}=32.9(\lambda L_{5100}/10^{44}\mathrm{ergs}\mathrm{s}^1)^{0.7}`$ light days (Kaspi et al. 2000). The FWHM of the H$`\beta `$ line is rather uncertain given the contamination from the strong stellar spectrum. The FWHM of Mg II is generally consistent with that of H$`\beta `$, certainly so within the uncertainty of this sort of black hole estimate, so we adopt FWHM $`(H\beta )`$ = 5500 km s-1(Brotherton et al. 1999a). This yields MBH = $`2.7\times 10^8`$ M, and MBH/MBulge = 0.04%, which is within the sample range of local AGN but smaller than the average $`<`$MBH/M$`{}_{bulge}{}^{}>`$ 0.1% (Merritt & Ferrarese 2001). We estimate the bolometric luminosity as $`L_{bol}=9\lambda L_\lambda `$, where $`\lambda L_\lambda `$ is evaluated at 5500 Å (Kaspi et al. 2000), obtaining $`5.7\times 10^{11}L_{}`$, consistent with our previous value (Brotherton et al. 1999a). For our black hole mass estimate, the Eddington luminosity is $`1.0\times 10^{13}L_{}`$, and the Eddington fraction $`L/L_{Edd}=0.05`$. These values are all consistent with those of other QSOs and Seyfert galaxies so far studied, so nothing stands out as particularly extreme other than the starburst.

While the extremely massive starburst of UN J1025$``$0040 might suggest that it is an isolated case of limited applicability, it actually represents the prototype of an entire class of post-starburst QSOs. The 2dF QSO Redshift Survey now contains close to 23000 QSOs at $`z<3`$ and $`B<21`$ (e.g., Croom et al. 2001). Some $``$3% of $`z<0.75`$ quasars (94/2939) show evidence of Balmer jumps and high-order Balmer absorption lines at the 3 $`\sigma `$ level (Croom et al. in prep.). Determining the starburst ages and masses in conjunction with the AGN properties (e.g., black hole mass, Eddington fraction as above) for a large sample may lead to a consistent picture for the relationship between AGNs and the evolution of the galaxies that host them, something that we cannot yet do given only the singular case of UN J1025$``$0040. For instance, we might expect to find that the starburst age correlates with MBH/Mbulge since older quasars will have had more time to accrete the available fuel, and inversely correlates with $`L/L_{Edd}`$ since that fuel must eventually be exhausted.

Quantitatively developing the evolutionary scenario will probably still be insufficient to explain all vectors leading to nuclear activity and their duty cycles. Besides the 97% of 2dF QSOs that do not show clear spectroscopic evidence for (extremely) massive starbursts and which may not represent evolved ULIRGs, there are binary quasars. There are around a dozen binary quasars currently known, and the indications are that the nuclei are triggered early in a galactic merger process by tidal forces when the galaxies are some 50-100 kpc apart (Mortlock et al. 1999). The post-starburst quasar that is FIRST J1643+3156B, the fainter member of a binary quasar system (Brotherton et al. 1999b), would fit right into the post-starburst quasar sample of the 2dF. High-quality follow-up imaging of that sample may reveal binary quasars like FIRST J1643+3156 in addition to evolved starburst mergers like UN J1025$``$0040. Starburst age would therefore appear to be a superior clock of fuel flow compared to merger age, in characterizing nuclear activity, and should be measured when possible.

Only a decade ago, debate focused upon whether the majority of quasars had “fuzz” around them or not. HST and ground-based near-IR imaging then opened up the study of quasar host galaxies in general. Today we are starting to characterize their stellar populations and may be on the verge of reconstructing their evolutionary histories. In the future, studies of the evolution of galaxies and AGN may be impossible to separate.

## 6 Summary

High-resolution images using WFPC2 on the Hubble Space Telescope show that UN J1025$``$0040 has a disturbed, compact nucleus, but no easily identifiable morphological structures. The photometry is consistent with that of Canalizo et al. (2000) indicating the AGN has not varied significantly between the time the ground-based data was taken and the HST images. The F450W and F814W images can be decomposed into a point source (the AGN) and residual, somewhat extended emission (starlight) concentrated within the central several hundred parsecs. The ratio of AGN to starlight in the red F814W image is consistent with the previous spectral analysis of Brotherton et al. (1999a), but the ratio in the blue F450W image is much smaller than expected, indicating the presence of a young stellar population. The currently available data constrain the young stars to be anywhere from a $``$40-Myr-old population with 10% of the mass of the dominant 400 Myr starburst, to a currently active starburst having 1% of the mass – these would still constitute massive starbursts in their own right. The older starburst dominates the light at red wavelengths, and the young starburst dominates the light at UV wavelengths. This young starburst may represent the trigger of the current AGN activity.

We thank Tod Lauer for providing his image-reduction software. Matthew Grabelsky acknowledges support from the NOAO/KPNO Research Experiences for Undergraduates program funded by the National Science Foundation. Support for proposal GO-08703 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. A.V.F. acknowledges the Guggenheim Foundation for a Fellowship. A portion of this work has been performed under the auspices of the U.S. Department of Energy, National Nuclear Security Administration by the University of California, Lawrence Livermore National Laboratory under Contract W-7405-ENG-48.