text:
# 1 Introduction

## 1 Introduction

Quantumchromodynamics (QCD) at finite temperature and chemical potential raises many interesting questions which may be relevant for the physics of the early universe, dense stars or relativistic heavy ion collisions. Since lattice calculations are not yet feasible at $`\mu 0`$, the QCD phase diagram is still to a large extent speculative, although there is little doubt that the phase structure is quite rich . In view of this situation, it may be of interest to go back to some exactly soluble field theoretic models and try to understand their phase diagram, in particular the structure of dense, baryonic matter (for a recent review, see Ref. ). The simplest such model is perhaps the Gross-Neveu model , $`N`$ species of massless fermions with a quartic self-interaction,

$$=\overline{\psi }\mathrm{i}\gamma ^\mu _\mu \psi +\frac{1}{2}g^2\left[(\overline{\psi }\psi )^2+\lambda (\overline{\psi }\mathrm{i}\gamma _5\psi )^2\right],\lambda =0,1.$$
(1)
(Flavor indices are suppressed, i.e., $`\overline{\psi }\psi =_{k=1}^N\overline{\psi }_k\psi _k`$ etc.) This model can be solved analytically in the large $`N`$ limit and shares some properties with QCD, notably asymptotic freedom. It exists in two variants with either discrete ($`\lambda =0`$) or continuous ($`\lambda =1`$) chiral symmetry. For the sake of simplicity, we shall refer to the first model as Gross-Neveu (GN) model, to the second one as 2-dimensional Nambu–Jona-Lasinio (NJL2) model . Aside from the spectrum which contains massive, unconfined fermions as well as composite mesons and baryons, the phase diagram of these models has also been studied some time ago . It is quite non-trivial, including a tricritical point separating lines of first and second order transitions in striking analogy with recent conjectures about real QCD  . Here, we shall be mainly interested in the behavior of these models at finite density but zero temperature. There, either the discrete or the continuous chiral symmetry which are spontaneously broken in the vacuum are believed to get restored in a first order phase transition at the critical density $`\rho _{\mathrm{crit}}=m/(\pi \sqrt{2})`$, with $`m`$ the dynamical fermion mass in the vacuum. Below this density, matter should consist of droplets of chirally restored phase containing extra fermions surrounded by regions of chirally broken vacuum, much like in the bag model. Above $`\rho _{\mathrm{crit}}`$ the fermions become massless and free.

This scenario which dates back to the mid 80’s was recently challenged in the case of the NJL2 model by one of us , and we extend our criticism to the GN model in the present work. Our main objection is very simple indeed: The conventional wisdom ignores the fact that these models possess baryons, either massless (NJL2 model) or massive ones (GN model) (see  and references therein). In the low density limit, one would expect that baryonic matter will form widely separated baryons, its energy density satisfying

$$\frac{}{\rho }|_{\rho =0}=M_B.$$
(2)
This relation is violated by the standard solutions to the GN or NJL2 models. The reason can be traced back to the tacit assumption of unbroken translational invariance inside and outside the droplets. In the case of the NJL2 model, the energetically favored structure of baryonic matter was shown to be a chiral crystal : Baryons are topological excitations of the pion field, analogous to the Skyrme model in higher dimensions , and matter is a kink crystal, the low dimensional analogue of a Skyrme crystal . This result can be shown analytically with the help of bosonization techniques; Eq. (2) is then indeed satisfied for the true ground state with $`M_B=0`$. Strictly speaking, no restoration of chiral symmetry takes place as a function of density, although a rapidly fluctuating condensate may be considered as equivalent to the chirally symmetric phase for most practical purposes. The same argument will apply for the crystalline configuration of the GN model constructed in this work.

The problem in the GN model is more involved, since the bosonization used in Ref.  can only be applied to the light “pion”. Therefore, we cannot present a solution to the finite density GN model (possessing only the massive “sigma” meson) which would be nearly as complete as for the NJL2 model. We have performed a variational calculation in order to find out whether the standard Fermi gas ground state of baryonic matter in the GN model is unstable against crystallization.

We finally remind the reader that quite generally, breaking of continuous symmetries and Goldstone bosons are not allowed in 1+1 dimensions . As is well known, the large $`N`$ limit suppresses the fluctuations which would otherwise restore the symmetry. This makes it possible to study such effects in 1+1 dimensions as well, provided one works only in leading order in $`1/N`$, i.e., takes the limit $`N\mathrm{}`$ at face value. This remark applies both to the continuous chiral symmetry of models like the NJL2 or ’t Hooft model  where massless, pion-like excitations appear, and to translational invariance in the GN model, our main theme here.

## 2 Construction of a trial wave function

In the large $`N`$ limit, the relativistic Hartree-Fock approximation for the ground state becomes exact. It corresponds to a variational calculation where the trial wave functions are all possible Slater determinants; the Hartree-Fock method determines the best fermion single particle orbits. Since we cannot solve this problem fully at this point, we search for the best Slater determinant in a restricted subspace only, generating the single particle orbits from a certain one-parameter family of periodic scalar potentials. These scalar potentials need not be self-consistent since they are only used to restrict the subspace of Slater determinants. By evaluating the expectation value of the exact Hamiltonian $`H`$ in this trial state and minimizing it with respect to this one parameter, we obtain an upper bound for the ground state energy. If this upper bound is below the energy of the translationally invariant Hartree-Fock solution, we can in this way establish instability of the system with respect to breakdown of translational invariance. In the present section, we outline the construction of our trial wave function. It is motivated primarily by our expectation that matter at low density will tend to form isolated baryons. Another guiding principle is the desire to gain as much analytical insight as possible.

Our starting point is the scalar potential for a single baryon in the Gross-Neveu model. Using units where $`m=1`$ ($`m`$ is the physical fermion mass), it is given by

$$S(x)=1+y(\mathrm{tanh}\xi _{}\mathrm{tanh}\xi _+),\xi _\pm =yx\pm \frac{1}{2}\mathrm{artanh}y.$$
(3)
The mean-field Dirac Hamiltonian with this potential,

$$h_D=\gamma _5\frac{1}{\mathrm{i}}\frac{}{x}+\gamma ^0S(x),$$
(4)
can be diagonalized analytically. The parameter $`y[0,1]`$ is related to the occupation of the valence level. If $`y\sqrt{3}/2`$, $`S(x)`$ crosses zero at the points $`x=\pm R/2`$ with

$$R=\frac{1}{2y}\mathrm{ln}\left(\frac{2y^21+yw}{2y^21yw}\right),w=\sqrt{4y^23}.$$
(5)
Furthermore, we note that as a result of spontaneous breakdown of the discrete $`\gamma _5`$ symmetry, $`S(x)`$ and $`S(x)`$ are equivalent scalar baryon potentials built on the two degenerate vacua. The corresponding single particle solutions of the Dirac equation can be related by $`\psi ^{}(x)=\mathrm{i}\gamma _5\psi (x)`$. This enables us to construct a family of smooth, periodic potentials for which the Dirac equation is again exactly solvable: Choose $`S(x)`$ from Eq. (3) in the interval $`[R/2,R/2]`$ and extend it to the whole axis via $`S(x+R)=S(x)`$. By pasting together pieces of the single baryon potential with alternating signs in this manner, we obtain a continuous and differentiable scalar potential whose period of $`2R`$ is controlled by the parameter $`y`$. The shape of the potential changes with the period in a prescribed way, see Fig. 1. In the limit $`y1`$, the potential corresponds to a succession of equidistant, alternating kinks and anti-kinks at a distance growing logarithmically with $`(1y)`$. In the opposite limit $`y\sqrt{3}/2`$, $`S`$ may be interpreted as arising from strongly overlapping baryons with an almost complete cancellation of their scalar potentials. Only some rapid residual oscillation around zero is left which goes over into the value zero of the chirally restored phase at the endpoint $`y=\sqrt{3}/2`$. The resulting scalar potential is not only periodic with period $`2R`$, but also invariant under translations by half a period ($`R`$) combined with a discrete $`\gamma _5`$-transformation. In this sense, it may be viewed as discrete version of the “chiral spiral” of Ref. . What we have set up here is evidently a kink-antikink crystal, as opposed to the kink crystal characteristic for the ground state of models with continuous chiral symmetry.

The Dirac equation for a single baryon can be solved in closed form . Since we shall be mostly interested in negative energy continuum states in what follows, we recall the well-known results for these particular wavefunctions only. In the interval $`[R/2,R/2]`$, one finds (we use a representation where $`\gamma ^0=\sigma _1,\gamma ^1=\mathrm{i}\sigma _3,\gamma _5=\gamma ^0\gamma ^1=\sigma _2`$),

$$\psi _k(x)=\left(\begin{array}{c}(\mathrm{i}k1)(\mathrm{i}ky\mathrm{tanh}\xi _{})\\ E(k)(\mathrm{i}ky\mathrm{tanh}\xi _+)\end{array}\right)\mathrm{e}^{\mathrm{i}kx},E(k)=\sqrt{k^2+1}.$$
(6)
$`\psi _k`$ and $`\psi _k`$ are degenerate with eigenvalue $`E(k)`$. The general solution of the stationary Dirac equation to $`E(k)`$ in this interval is then

$$\psi (x)=A\psi _k(x)+B\psi _k(x).$$
(7)
Furthermore there exists a discrete valence level at energy $`\sqrt{1y^2}`$.

Let us now apply the textbook procedure from solid state physics in order to extend the solution to the whole $`x`$-axis and determine the spectrum of the Dirac Hamiltonian. In a slight modification of Bloch’s theorem, we require the following boundary condition,

$$\psi (R/2)=\mathrm{e}^{\mathrm{i}\delta }\mathrm{i}\gamma _5\psi (R/2).$$
(8)
This guarantees a smooth matching of the wave functions at the borders of different intervals. Inserting the ansatz (7) and denoting the upper (lower) components of $`\psi _k`$ by $`\phi _k`$ ($`\chi _k`$), Eq. (8) yields the generalized eigenvalue problem

$$Mv=\mathrm{e}^{\mathrm{i}\delta }Nv$$
(9)
with

$$M=\left(\begin{array}{cc}\phi _k(R/2)& \phi _k(R/2)\\ \chi _k(R/2)& \chi _k(R/2)\end{array}\right),N=\left(\begin{array}{cc}\hfill \chi _k(R/2)& \hfill \chi _k(R/2)\\ \hfill \phi _k(R/2)& \hfill \phi _k(R/2)\end{array}\right),v=\left(\begin{array}{c}A\\ B\end{array}\right).$$
(10)
Hence the factors $`\mathrm{e}^{\mathrm{i}\delta }`$ are the eigenvalues of the $`2\times 2`$-matrix $`N^1M`$. Evaluation of this matrix using Eq. (6) shows that it has unit determinant, whereas its trace is given by

$$\mathrm{tr}\left(N^1M\right)=2X_k,$$
(11)
with

$$X_k=\frac{E(k)}{k(k^2+y^2)}\left[(k^2+1y^2)\mathrm{sin}(kR)+kw\mathrm{cos}(kR)\right].$$
(12)
The characteristic equation then reads

$$(\mathrm{e}^{\mathrm{i}\delta })^22X_k\mathrm{e}^{\mathrm{i}\delta }+1=0.$$
(13)
For $`|X(k)|1`$, the matrix $`N^1M`$ is unitary, $`\delta =\mathrm{arccos}X_k`$ is real and the corresponding level lies inside an allowed band. For $`|X(k)|>1`$, $`\delta `$ is imaginary and we set $`\delta =\mathrm{i}\zeta `$. The eigenvalues of $`N^1M`$ are now $`\mathrm{e}^{\pm \zeta }`$ with wave functions blowing up in one direction — the corresponding state is forbidden and lies inside a band gap. The diagonalization of Eq. (9) also determines the ratio $`A/B`$ of the coefficients introduced in Eq. (7); one finds

$$\frac{A}{B}=\frac{1+\mathrm{i}k}{1y^2}\left(Y_k\pm \sqrt{Y_k^2(1y^2)^2/E(k)^2}\right),$$
(14)
with

$$Y_k=kw\mathrm{sin}(kR)(k^2+1y^2)\mathrm{cos}(kR).$$
(15)
The two solutions correspond to the two degenerate states with energy $`E(k)`$. This procedure yields both the spectrum of the Dirac Hamiltonian and the single particle wave functions up to the normalization for which we shall use a box normalization later on (cf. Sect. 3). Here, let us briefly illustrate the resulting band structure of the potential used to generate our trial wave function. In Fig. 2, we show an example of the single particle spectrum, comparing it to the one corresponding to the interacting vacuum (or, equivalently, free massive fermions with $`m=1`$). It exhibits the expected band structure. Since all eigenvalues still occur in pairs $`\pm E(k)`$, the Dirac sea inherits this band structure, a novel feature as compared to a non-relativistic solid. Fig. 3 displays the lowest few bands as function of the parameter $`y`$. The left-hand part of this figure (corresponding to small lattice spacing $`2R`$) can be understood perturbatively, since the scalar potential is weak (cf. Fig. 1). Using degenerate perturbation theory, one easily finds the following position and width of the gaps, expressed in terms of $`\eta =y\sqrt{3}/2`$,

$$E_{\mathrm{gap}}\frac{(2n+1)\pi }{3^{(1/4)}8\sqrt{\eta }}\left(1\frac{\eta }{2\sqrt{3}}\right),\mathrm{\Gamma }_{\mathrm{gap}}\frac{64\sqrt{3}\eta }{\pi ^3(2n+1)^3},n=0,1,2\mathrm{}$$
(16)
On the right-hand side of Fig. 3 (large lattice spacing), all gaps move into the mass gap and we recover the spectrum of the Callan-Coleman-Gross-Zee kink , including the valence state at zero energy. The lowest allowed band obviously develops out of the valence level of the single baryon (which moves to zero energy as $`y1`$), see the dashed curve in Fig. 3. Finally, in Fig. 4, we illustrate the band structure together with the corresponding scalar potential, for three values of $`y`$. Case a) features well separated baryons, weak tunneling and hence a narrow valence band. In case b), the band width is comparable to the depth of scalar potential. In case c), the valence band has become much wider than the amplitude of the scalar potential. Fermions can move freely through the whole crystal and baryons have lost their identity. By varying the parameter $`y`$, we let the system decide which option it prefers at a given density.

## 3 Computing the regularized ground state energy

In our variational calculation of the ground state, we proceed as follows: We take the periodic scalar potential from the preceding section and fill the single particle states until the desired average baryon density is reached. Next, we evaluate the expectation value of the exact Hamiltonian in this trial state and minimize with respect to $`y`$, our only variational parameter. At the minimum, there will be a unique relationship between lattice spacing $`2R`$ (or, equivalently, $`y`$) and baryon density. In the case of the single baryon, $`y`$ was linked to the occupation fraction $`\nu =n/N`$ of the valence level via $`y=\mathrm{sin}(\pi \nu /2)`$.

Evaluation of the ground state energy cannot be done naively for several reasons. First, there is a trivial quadratic divergence in the kinetic energy due to the Dirac sea. This will be removed by subtracting the ground state energy of a free Fermi gas of corresponding density. Then, there are logarithmic divergences both in the kinetic and in the potential energy which require the usual renormalization of the theory at zero density. Moreover, due to the gap structure of the Dirac sea, a direct calculation in the continuum would be technically difficult. We therefore enclose the system in a box of size $`L`$ (an integer multiple of the lattice spacing, $`L=2RK`$) and require antiperiodic boundary conditions, thereby discretizing the spectrum of $`h_D`$. The size of the box will be made large enough so that our results do not depend on the IR-regulator $`L`$. The phase factor $`\mathrm{e}^{\mathrm{i}\delta }`$ in Eq. (8) then has to satisfy

$$(\mathrm{e}^{2\mathrm{i}\delta })^K=1,\delta =p_nR+\frac{\pi }{2},p_n=\frac{\pi }{L}(2n+1)$$
(17)
with $`n`$ an integer. Solving the characteristic equation (13) yields

$$\mathrm{e}^{\mathrm{i}\delta }=X_k\pm \mathrm{i}\sqrt{1X_k^2}$$
(18)
or, equivalently,

$$\mathrm{sin}(p_nR)=X_k.$$
(19)
A solution of Eq. (19) for $`k`$ then yields the single particle energies $`\pm E(k)`$ corresponding to the discrete $`p_n`$.

To reach a certain baryon density, one has to fill all negative energy states plus a number of positive energy states. “Antimatter” on the other hand is obtained by leaving a number of negative energy levels unoccupied. Since the energy density is identical in both cases, we have found it more convenient to lower the level of the Dirac sea rather than raise it. This is the reason why — as mentioned above — we need to consider only negative energy single particle states below the “Fermi surface”. We denote the single particle wave functions by $`\psi _{n,i}`$, where $`n`$ refers to the energy and $`i=1,2`$ labels the two degenerate solutions.

Taking the expectation value of $`H`$ in our trial state and invoking the large $`N`$ limit yields the following expression for the ground state energy per flavor,

$$E=T+V$$
(20)
with

$$T=\underset{n}{\overset{occ}{}}\underset{i=1}{\overset{2}{}}2K_{R/2}^{R/2}dx\frac{1}{𝒩_{n,i}}\psi _{n,i}^{}\gamma _5\frac{1}{\mathrm{i}}\frac{}{x}\psi _{n,i},$$
(21)
$$V=Ng^2K_{R/2}^{R/2}dx\overline{\psi }\psi ^2,\overline{\psi }\psi =\underset{n}{\overset{occ}{}}\underset{i=1}{\overset{2}{}}\frac{1}{𝒩_{n,i}}\overline{\psi }_{n,i}\psi _{n,i}.$$
(22)
The normalization factors are

$$𝒩_{n,i}=2K_{R/2}^{R/2}dx\psi _{n,i}^{}\psi _{n,i}.$$
(23)
The kinetic energy has both quadratic and logarithmic divergences. Let us isolate the divergences by adding and subtracting the asymptotic form ($`n\mathrm{}`$) of the terms in Eq. (21),

$$T=T_{\mathrm{conv}}+T_{\mathrm{div}},$$
(24)
with

$`T_{\mathrm{conv}}`$ $`=`$ $`{\displaystyle \underset{n}{\overset{occ}{}}}\left\{{\displaystyle \underset{i=1}{\overset{2}{}}}2K{\displaystyle _{R/2}^{R/2}}dx{\displaystyle \frac{1}{𝒩_{n,i}}}\psi _{n,i}^{}\gamma _5{\displaystyle \frac{1}{\mathrm{i}}}{\displaystyle \frac{}{x}}\psi _{n,i}+2p_n{\displaystyle \frac{1}{p_n}}\left(1{\displaystyle \frac{2w}{R}}\right)\right\},`$
$`T_{\mathrm{div}}`$ $`=`$ $`{\displaystyle \underset{n}{\overset{occ}{}}}\left\{2p_n+{\displaystyle \frac{1}{p_n}}\left(1{\displaystyle \frac{2w}{R}}\right)\right\}.`$ (25)
Here, $`p_n`$ are the momenta of the free fermions introduced in Eq. (17). One can derive this asymptotic form using the known single particle wave functions of Sect. 2. The potential energy has a logarithmic divergence which can be exposed as follows. First, split the condensate into regular and singular pieces,

$$\overline{\psi }\psi =\overline{\psi }\psi _{\mathrm{conv}}+\overline{\psi }\psi _{\mathrm{div}},$$
(26)
with

$`\overline{\psi }\psi _{\mathrm{conv}}`$ $`=`$ $`{\displaystyle \underset{n}{\overset{occ}{}}}{\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{1}{𝒩_{n,i}}}\overline{\psi }_{n,i}\psi _{n,i}+{\displaystyle \frac{1}{L}}{\displaystyle \underset{n}{\overset{occ}{}}}\left({\displaystyle \frac{2S}{p_n}}\right),`$
$`\overline{\psi }\psi _{\mathrm{div}}`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{n}{\overset{occ}{}}}\left({\displaystyle \frac{2S}{p_n}}\right).`$ (27)
Once again, the explicit form of the single particle wave functions has been used to derive Eq. (27), where $`S`$ is the scalar potential described above. As usual, the gap equation at zero density is used to renormalize the coupling constant,

$$\frac{1}{Ng^2}=\frac{2}{L}\underset{n}{\overset{vac}{}}\frac{1}{E(p_n)}.$$
(28)
Here, the summation over $`n`$ is according to the vacuum occupation, as opposed to the occupation for baryonic matter which enters the sums in Eqs. (21,22). The UV-cutoff $`\mathrm{\Lambda }`$ has to be the same in both sums. We then insert the decomposition (26) into Eq. (22). Since $`Ng^2`$ vanishes like $`(\mathrm{ln}\mathrm{\Lambda })^1`$, we may drop the $`\overline{\psi }\psi _{\mathrm{conv}}^2`$ term and have to keep only the two remaining terms,

$$V=V_1+V_2$$
(29)
with

$`V_1`$ $`=`$ $`Ng^2K{\displaystyle _{R/2}^{R/2}}dx\mathrm{\hspace{0.17em}2}\overline{\psi }\psi _{\mathrm{conv}}\overline{\psi }\psi _{\mathrm{div}},`$
$`V_2`$ $`=`$ $`Ng^2K{\displaystyle _{R/2}^{R/2}}dx\overline{\psi }\psi _{\mathrm{div}}^2.`$ (30)
Using Eqs. (27) and (28) and performing the limit $`\mathrm{\Lambda }\mathrm{}`$, $`V_1`$ reduces to

$$V_1=2K_{R/2}^{R/2}dx\overline{\psi }\psi _{\mathrm{conv}}S.$$
(31)
$`V_2`$ is given by

$$V_2=\frac{\left(_n^{occ}\frac{1}{p_n}\right)^2}{\left(_n^{vac}\frac{1}{E(p_n)}\right)}\frac{1}{R}_{R/2}^{R/2}dxS^2.$$
(32)
The integral can be performed in closed form with the compact result

$$_{R/2}^{R/2}dxS^2=(R2w).$$
(33)
Assuming the cutoff $`\mathrm{\Lambda }`$ to be large then yields

$$V_2=\left\{\underset{n}{\overset{occ}{}}\frac{1}{p_n}\left(\underset{n}{\overset{vac}{}}\frac{1}{E(p_n)}\underset{n}{\overset{occ}{}}\frac{1}{p_n}\right)\right\}\left(1\frac{2w}{R}\right).$$
(34)
If we add up $`T_{\mathrm{div}}`$, Eq. (25) and $`V_2`$, Eq. (34), and eliminate the trivial quadratic divergence by subtracting the free Fermi gas kinetic energy at the same density, we obtain the UV finite result

$$T_{\mathrm{div}}+V_2+\underset{n}{\overset{occ}{}}2p_n=\left(1\frac{2w}{R}\right)\left(\underset{n}{\overset{vac}{}}\frac{1}{E(p_n)}\underset{n}{\overset{occ}{}}\frac{1}{p_n}\right).$$
(35)
After some rearrangement, the total regularized energy can be cast into the form

$`E_{\mathrm{reg}}`$ $`=`$ $`{\displaystyle \underset{n}{\overset{occ}{}}}\left\{{\displaystyle \underset{i=1}{\overset{2}{}}}2K{\displaystyle _{R/2}^{R/2}}dx{\displaystyle \frac{1}{𝒩_{n,i}}}\psi _{n,i}^{}\gamma _5{\displaystyle \frac{1}{\mathrm{i}}}{\displaystyle \frac{}{x}}\psi _{n,i}+2p_n{\displaystyle \frac{1}{p_n}}\left(1{\displaystyle \frac{2w}{R}}\right)\right\}`$ (36)
$`+2K{\displaystyle _{R/2}^{R/2}}dxS\left\{{\displaystyle \underset{n}{\overset{occ}{}}}{\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{1}{𝒩_{n,i}}}\overline{\psi }_{n,i}\psi _{n,i}+{\displaystyle \frac{S}{L}}\left({\displaystyle \underset{n}{\overset{vac}{}}}{\displaystyle \frac{1}{E(p_n)}}+{\displaystyle \underset{n}{\overset{occ}{}}}{\displaystyle \frac{1}{p_n}}\right)\right\}.`$
Invoking the fact that the $`\psi _{n,i}`$ are eigenfunctions of the Dirac Hamiltonian (4) with eigenvalues $`E(k_n)`$ as well as the result (33), this expression can further be simplified to

$$E_{\mathrm{reg}}=2\underset{n}{\overset{occ}{}}\left(E(k_n)p_n\right)+\left(1\frac{2w}{R}\right)\underset{n}{\overset{vac}{}}\frac{1}{E(p_n)}.$$
(37)
The $`k_n`$ are the single particle momenta to the corresponding discrete $`p_n`$, obtained via the transcendental equation (19).

Remarkably, the total energy has been expressed entirely in terms of the single particle energies and known functions of $`y`$. That the sum in $`E_{\mathrm{reg}}`$ is UV finite can now be confirmed by means of the asymptotic behavior

$$k_np_n\frac{w}{Rp_n}(n\mathrm{}),$$
(38)
which in turn follows from the transcendental equation (19). Summarizing, the calculation of the variational energy has been reduced to solving a transcendental equation for $`k_n`$, Eq. (19), and performing a convergent sum, Eq. (37). The baryon density is adjusted via the occupation of the single particle states. There is only one variational parameter $`y`$ — it controls the spatial period of the scalar potential underlying our trial wave function.

## 4 Results of the variational calculation

We have performed calculations along the lines outlined above and evaluated the ground state energy of baryonic matter with our variational ansatz, the kink-antikink crystal. The size of the box $`L`$ was increased until the results were stable which required values of $`L3001600`$ (in units where $`m=1`$) increasing with density. The sum to be performed is UV finite so that no cutoff is required. To get precise results, we had to take into account single particle orbits with energies up to $`|E|50100`$ typically.

We first summarize our findings concerning the variational parameter $`y`$. The result is extremely simple and can be stated as follows: At all densities considered, the best choice of $`y`$ was the one which yielded a fully occupied valence band in baryonic matter (or empty in “antimatter”), cf. the shading in Fig. 4. With hindsight, this is rather plausible: Around each gap, the levels are pushed apart symmetrically due to the interaction (in first order perturbation theory). When summing over occupied states the main effect cancels out. Only at the Fermi surface the filling is asymmetric and the system can take advantage of the level shifts to lower its energy. This picture was confirmed by studying in detail where the bulk of the energy difference came from. It was indeed the region in the vicinity of the Fermi surface. The optimal value of $`y`$ translates into the following relation between baryon density and baryon radius (or lattice spacing),

$$\rho _B=\frac{1}{2R}.$$
(39)
By this relation, the different values for $`y`$ in Fig. 4 correspond to baryon densities a) $`\rho =0.12`$, b) $`\rho =0.3`$, and c) $`\rho =1.0`$. Each potential well may be associated with one baryon, although for high densities fermions are not confined to one well, cf. Fig. 4.

With this optimal choice of $`y`$ we get a unique answer for the best ground state energy as a function of density. It is shown in Fig. 5, together with the energies of the free Fermi gas and the energy for the standard “mixed phase” solution. We recall that in the latter phase, the extra fermions are put into droplets of chirally restored phase, like in the bag model, whereas the vacuum in between the baryons remains in the broken phase. The standard first order phase transition occurs at $`p_f=1/\sqrt{2}`$, or $`\rho 0.225`$, as indicated by the dot in Fig. 5. Our result is always below these other two curves. At small densities, it perfectly matches the expected slope, Eq. (2), with $`M_B=2/\pi `$, the mass of a single baryon with fully occupied valence level. This is no surprise since our variational ansatz was designed to have this limit correctly built in, but it is a good test of the whole calculation, including the regularization and renormalization procedure. In our ansatz, we had to assume that the distance between two baryons equals the baryon diameter in order to generate a smooth periodic potential by gluing together pieces of the single baryon potentials. This is certainly the right thing to do at very low densities. There, the extra fermions are concentrated in the kinks and anti-kinks. The interaction between kinks and anti-kinks is repulsive, therefore they will tend to spread equally; there is no distinction between left and right neighbours. It is more surprising that such an ansatz is also capable of lowering the energy at higher densities, even at densities where the standard picture predicts a transition to the chirally restored phase. Thereby, this standard first-order phase transition in the GN model can be ruled out, because our variational upper bound for the ground state energy excludes the chirally symmetric solution. In a sense, the mixed phase picture may be regarded as a crude attempt to describe baryons. It misses some important binding effects which survive the large $`N`$ limit , as manifested through the existence of single baryons in the spectrum of the GN model.

We have pushed the calculation up to higher densities without any indication that the curves shown in Fig. 5 will cross each other, see Fig. 6. The energy difference becomes tiny and the volume one needs to get the same accuracy becomes huge, so that it gets increasingly difficult to get precise numerical results. Our trial state definitely becomes worse with increasing density, as can be judged by comparing the chiral condensate $`\overline{\psi }\psi `$ computed from our wave functions with the input scalar potential. Hence the true ground state must be lower in energy at high densities than suggested by our variational calculation.

Finally, in Fig. 7, we illustrate the spatial dependence of the baryon density for the three values of $`y`$ already used in Fig. 4. Unlike in the NJL2 model where the baryon density was constant, here it also reflects the periodicity of the scalar potential.

Needless to say, it would be of quite some theoretical interest to find the truly self-consistent periodic Hartree-Fock potential for the GN model at finite density. This may be the prerequisite for revising the phase diagram of the GN model in the whole $`(T,\mu )`$ plane.


format: text

text:
# Quantum Phase Transition in Coupled Spin Ladders

\[

## Abstract

The ground state of an array of coupled, spin-half, antiferromagnetic ladders is studied using spin-wave theory, exact diagonalization (up to 36 sites) and quantum Monte Carlo techniques (up to 256 sites). Our results clearly indicate the occurrence of a zero-temperature phase transition between a Néel ordered and a non-magnetic phase at a finite value of the inter-ladder coupling ($`\alpha _c0.3`$). This transition is marked by remarkable changes in the structure of the excitation spectrum.

\]

Ground-state (GS) correlations in a quantum antiferromagnet are closely related to the nature of low-energy excitations. This is deeply connected to the mechanism of spontaneously symmetry breaking and turns out clearly whenever competing GS’s give rise to a quantum phase transition, leading in general to remarkable changes of the excitation spectrum at the critical point .

A simple model system experiencing such changes in the structure of low-energy excitations in correspondence of a quantum phase transition is the spin-half Heisenberg antiferromagnet on a two-dimensional array of coupled ladders:

$$\widehat{}=J\underset{r}{}\widehat{𝐒}_r\widehat{𝐒}_{r+\widehat{x}}+\underset{r}{}J_r\widehat{𝐒}_r\widehat{𝐒}_{r+\widehat{y}},$$
(1)
where $`\widehat{𝐒}_r=(\widehat{S}_r^x,\widehat{S}_r^y,\widehat{S}_r^z)`$ are $`s=1/2`$ operators on the sites $`r=(r_x,r_y)`$ of a $`L\times L`$ lattice with periodic boundary conditions; $`\widehat{x}=(1,0)`$, $`\widehat{y}=(0,1)`$, $`J_r=J`$ or $`J_r=\alpha J`$ ($`J>0`$), depending on the parity of $`r_y`$, and $`\alpha `$ is the inter-ladder coupling (see Fig. 1). Such Hamiltonian interpolates between the Heisenberg model on the square lattice ($`\alpha =1`$) and a system of $`L/2`$ decoupled two-leg ladders ($`\alpha =0`$). In the square lattice limit, the GS has Néel long-range order, with a gapless excitation spectrum and a sizable value of the antiferromagnetic order parameter . The two-leg ladder, instead, has a finite triplet gap in the thermodynamic limit and no long-range antiferromagnetic order . As a result a quantum critical point between a gapless, magnetically ordered phase and a non-magnetic GS of purely quantum mechanical nature is expected at a critical value of the inter-ladder coupling, $`\alpha _c`$ .

Besides its intrinsic theoretical interest, this model has been recently studied due to the discovery of several compounds, such as $`\mathrm{SrCu}_2\mathrm{O}_3`$ and $`(\mathrm{VO})_2\mathrm{P}_2\mathrm{O}_7`$, displaying clear signatures of a finite gap ($`J/2`$) in the excitation spectrum related to the underlying ladder structure . In addition, the coupled-ladder Hamiltonian (1) has been also considered as a simplified model for the magnetism of the striped phase in $`\mathrm{CuO}_2`$ planes of hole-doped $`highT_c`$ copper-oxides .

In this paper, by means of spin-wave (SW) theory, exact diagonalization by the Lanczos method and zero-temperature quantum Monte Carlo techniques, including variational and Green function Monte Carlo (GFMC) , we focus on the GS and the low-energy excitations of the coupled spin-ladder model (1). In particular, by means of a systematic size-scaling analysis, we will show how the structure of low-energy excitation spectrum provides very clear indications of the changes of the GS state correlations occurring at a critical point, thus allowing us to put on firmer grounds the existence of a quantum phase transition at $`\alpha _c0.3`$.

The simplest approach to the study the effects of zero-point fluctuations on the GS correlations of a quantum magnet is SW theory. This has turned out to be a very reliable approximation of the GS of spin-half systems whenever it has long-range antiferromagnetic order even in presence of strong quantum fluctuations . By the use of the standard Holstein-Primakoff representation of the spin operators we can compute the fluctuations over the classical solution at the leading order in $`1/s`$. In contrast to the Heisenberg model on the square lattice, the reduced translation symmetry along the $`y`$ direction, implies the existence of two inequivalent branches of SW excitations, and the energy of such SW modes, $`\omega _k^\pm `$, reads

$$\omega _k^\pm =[2D^2(\gamma _k^2+\gamma _{\overline{k}}^2)2\delta _k^2\pm F_k]^{1/2}/\sqrt{2},$$
(2)
with $`F_k=[(\gamma _k^2\gamma _{\overline{k}}^2)^2+4\delta _k^2(\gamma _k+\gamma _{\overline{k}})^2]^{1/2}`$, $`\gamma _k=(\beta \mathrm{cos}k_x+\mathrm{cos}k_y)/2`$, $`\delta _k=(1\alpha )/2(1+\alpha )\mathrm{sin}k_y`$, $`\overline{k}=k+(0,\pi )`$, $`D=(1+\beta )/2`$, and $`\beta =2/(1+\alpha )`$. The dispersion relation of the two SW modes is plotted in Fig. 2. Notice that for $`\alpha <1`$ a finite gap develops between the optical and the acoustical branch of the SW dispersion along the $`y`$-direction due to the reduced translation symmetry. The acoustical branch remains always gapless at $`k=(0,0)`$ and $`k=(\pi ,0)`$, corresponding to the Goldstone modes (in a reduced-zone scheme) associated with the SU(2) symmetry-breaking assumption.

The expansion of the staggered magnetization at the first order in $`1/s`$ indicates the stability of the Néel order up to very small values of the inter-ladder coupling ($`\alpha _c0.01`$, for $`s=1/2`$). Within the same approximation the GS energy per site reads

$$\frac{E_0^{SW}}{L^2}=\stackrel{~}{J}s(s+1)(1+\beta )+\frac{2\stackrel{~}{J}s}{L^2}\underset{k}{\overset{RBZ}{}}(\omega _k^++\omega _k^{}),$$
(3)
where $`\stackrel{~}{J}=(1+\alpha )J/2`$ and the summation is over the $`k`$-vectors belonging to the sector of the Brillouin zone with $`\pi /2<k_y\pi /2`$. The SW energy in the thermodynamic limit as a function of the inter-ladder coupling $`\alpha `$ is plotted in Fig. 3-a. For spin models the exact value of the GS energy can be calculated numerically for small clusters with the Lanczos technique and for rather larger sizes, in absence of frustration, with the GFMC method . The comparison with the extrapolation to the bulk limit of the GFMC results (up to $`L16`$, see below for details) indicates that the latter analytical approach provides accurate estimates of the GS energy only in the regime of large inter-ladder coupling ($`\alpha 0.6`$) thus suggesting that the SW approach underestimates the actual effect of quantum fluctuations and, therefore, the critical value of the inter-ladder coupling.

By means of the SW analysis it is also possible to derive a variational wave function providing a good representation of the GS and low-lying excited states at least for $`\alpha 1`$, when long-range antiferromagnetic correlations are expected. This wave function, which is also easily computable when used for importance sampling in a GFMC calculation to reduce the numerical effort , can be obtained starting from a Néel ordered state and including Gaussian fluctuations by means of a Jastrow factor :

$$|AF=𝒫_S\underset{x}{}S_M(x)\mathrm{exp}\left[\frac{1}{2}\underset{r,r^{}}{}v(rr^{})S_r^zS_r^{}^z\right]|x.$$
(4)
Here $`|x`$ is an Ising spin configuration specified by assigning the value of $`S_r^z`$ for each site, $`𝒫_S`$ is the projector onto the subspace with $`S_{tot}^z=_rS_r^z=S`$, and $`S_M(x)=(1)^{N_{}(x)}`$ is the Marshall sign, reproducing exactly the phases of the GS . This depends only on the number of up spins on one of the two sublattices, $`N_{}(x)`$, so that $`|N=_xS_M(x)|x`$ represents the classical Néel state . For the two-body Jastrow potential, the simple form, based on the consistency with linear SW theory in the square lattice case, $`v(r)=(\eta /L^2)_{k0}e^{ikr}v_k`$, with $`v_k=1\sqrt{(1+\mathrm{\Gamma }_k)/(1\mathrm{\Gamma }_k)}`$, $`\mathrm{\Gamma }_k=(\mathrm{cos}k_x+\mathrm{cos}k_y)/2`$, and $`\eta `$ variational parameter, provides good variational estimates only for values of $`\alpha `$ very close to 1 (see Fig. 3-b), as it is expected since the modulation of the exchange interaction weakens the antiferromagnetic ordering and enhances spin fluctuations.

Fairly better variational results can be obtained by generalizing the latter wave function to the case of coupled spin ladders using the results of the SW analysis presented above. Indeed, the resulting wave function, whose explicit expression is too cumbersome to be reported here, provides very good variational estimates up to remarkably smaller values of $`\alpha 0.6`$ (Fig. 3-b). Decreasing further the inter-ladder coupling the accuracy of the SW variational wave function rapidly decreases. This clearly indicates that entering the small-$`\alpha `$ regime our gapless Néel ordered variational ansatz does not reproduce correctly the GS correlations.

Useful indications on the nature of the thermodynamic GS can be obtained numerically by studying the finite-size scaling of the GS energy, which is deeply connected to the nature of the excitation spectrum and therefore to possible thermodynamic broken symmetries . In fact, in presence of Néel long-range order, being the spectrum gapless and the magnon dispersion relation linear in the wave vector $`k`$, the leading finite-size correction to the GS energy per site, $`e_0(L)=E_0(L)/L^2`$, is $`𝒪(L^3)`$. Instead, in presence of a finite correlation length, a finite gap in the spin excitation spectrum and an exponential asymptotic dependence is expected. Here we have assumed a size dependence of the form

$$e_0(L)=e_0(\mathrm{})+c_0\mathrm{exp}(L/L_0)/L^2,$$
(5)
which has been employed in Ref.  for the two-leg ladder. As it is shown in Fig. 4, the size-scaling law predicted for a long-range ordered GS is fulfilled for values of $`\alpha `$ close enough to 1 while for small value of $`\alpha `$ a clear deviation from this behavior is observed and the exponential law (5) is instead satisfied. This provides another numerical evidence of the melting of the antiferromagnetic long-range order due to the transition to the ladder-like regime.

The occurrence of a quantum phase transition to a non-magnetic GS is clearly confirmed by the study of the spin gap shown in Fig. 5. This physical quantity can be measured straightforwardly with GFMC by performing two different simulations in the $`S_{tot}^z=0`$ and $`S_{tot}^z=1`$ subspaces. In contrast to finite-temperature algorithms, the calculation of the spin gap within GFMC does not involve any fitting procedure of low-temperature data and it is exact within statistical error. As shown in the left panel of Fig. 5, the triplet gap, $`\mathrm{\Delta }=E_1E_0`$, is a decreasing function of the inter-ladder coupling even if, on finite-sizes, this quantity is always non-zero for any value of $`\alpha `$. However, the size scaling of the spin gap, shown in the right panel of the same figure, indicates a clear deviation from the size dependence expected in presence of long-range Néel order ,

$$\mathrm{\Delta }_L=a/L^2+b/L^3,$$
(6)
and the opening of a finite gap in the thermodynamic excitation spectrum for $`\alpha 0.35`$. In this regime the extrapolation to the bulk limit of the finite-size data can be done using a law of the type:

$$\mathrm{\Delta }_L=\mathrm{\Delta }+a/L^2+b/L^4.$$
(7)
An estimate of the critical value of the inter-ladder coupling can be obtained by fitting the extrapolated values of the gap with the scaling law $`\mathrm{\Delta }=(\alpha _c\alpha )^\nu `$, with $`\nu 0.69`$, predicted for a quantum phase transition in (2+1) dimensions. This procedure gives the value $`\alpha _c0.32\pm 0.03`$, in agreement with previous numerical estimates obtained with finite-temperature algorithms  and the mean-field predictions of Ref.  but in contrast with the conclusions of bond-mean-field theory  indicating the vanishing of the gap for infinitesimal values $`\alpha `$. In contrast to frustrated systems like the $`J_1J_2`$ model, where the SW theory provides an accurate prediction of the transition to a non-magnetic phase due to competing interactions , in this case the SW result ($`\alpha _c0.01`$) grossly underestimate the exact one. This can be ascribed to the fact that in this model there are no competing GS’s at the classical level, the Néel state being stable up to $`\alpha =0`$, and the transition has therefore a pure quantum origin. Within the $`1/s`$ expansion, instead, the SW velocity remains finite up to $`\alpha =0`$ and the reduction of the staggered magnetization is only due to the crossover to the one-dimensional regime in which antiferromagnetic long-range order is unstable.

A major fingerprint of the dramatic changes in the structure of the excitation spectrum occurring in correspondence of the above quantum phase transition can be found in the finite-size behavior of the so called rigid rotator anomaly , $`\delta \chi ^1(S)=1/2\chi _S1/2\chi _L`$, with

$$\frac{1}{2\chi _S}=L^2\frac{E_SE_0}{S(S+1)},$$
(8)
where $`E_S`$ is the energy of the lowest excitation with spin $`S`$. In fact, in presence of long-range antiferromagnetic order, the low-lying excited states of spin $`S`$ are predicted to behave as the spectrum of a free quantum rotator, $`E_SE_0S(S+1)/L^2`$, as long as $`SL^2`$. In contrast, the expected behavior for a spin ladder is $`E_SE_0S`$, as it is easy to understand in a spin-liquid Resonating Valence Bond (RVB) picture. As a result, in the gapless phase the rigid rotator anomaly $`\delta \chi ^1`$ has to vanish identically in the thermodynamic limit while in the gapped regime $`\delta \chi ^1(S)L(L/(S+1)1)`$ diverges linearly with the volume. These features of the excitation spectrum clearly discriminate the two different zero-temperature phases of the present model, as it is evident from Fig. 6 displaying the very different behavior of $`\delta \chi ^1`$ below and above the critical point.

In summary, we have investigated the ground-state properties of an array of coupled spin-half ladders using spin-wave theory and a numerical analysis of the finite-size low-energy excitation spectrum. While the spin-wave theory turns out to be reliable only in the regime of weak inter-ladder coupling, our numerical results provide robust indications of the opening of a finite spin gap in the thermodynamic limit for $`\alpha 0.3`$, corresponding to a quantum phase transition between a gapless Néel ordered phase and a spin-liquid RVB ground state.

We are deeply indebted to Alberto Parola and Sandro Sorella for many valuable discussions and for their constant encouragement. Thanks also to Frederic Mila and Ingo Peschel for suggestions and useful discussions. This work was partially supported by MURST (COFIN00).


text:
# Testing the Supernova, Cepheid, and Early-type Galaxy Distance Scales

## 1. Early-type Galaxy Comparisons: SBF versus FP

The two most frequently applied early-type galaxy distance indicators are the fundamental plane (FP, and the related $`D_n`$-$`\sigma `$) and surface brightness fluctuations (SBF) methods. In a recent study (Blakeslee et al. 2001, 2002), we used $`V`$\- and $`I`$-band data from the ground-based SBF Survey (Tonry et al. 2001) to calculate FP photometric parameters for 170 galaxies with velocity dispersions available in the homogenized SMAC catalogue (Hudson et al. 2001). To our knowledge, this is the largest galaxy-by-galaxy comparison of different standard candle/rod distance methods to date. Fig. 1a shows the comparison.

Overall the distance agreement was good, but several low-luminosity, S0 galaxies had systematically low FP distances, probably due in part to younger ages and lower mass-to-light ratios, although aperture effects may also contribute. The SBF distances are tied to the Cepheids via measurements in spiral bulges, while the FP distances are tied to the Hubble flow via distant clusters; the Hubble constant that results from this comparison is $`H_0=68`$ km s$`^1`$Mpc<sup>-1</sup>. However, we also derived independent distances for these galaxies based on their velocities and the gravity field inferred from the redshift-space galaxy density; the resulting comparison with SBF yields $`H_0=74`$ (Fig. 1b), formally discordant at the 2$`\sigma `$ level with the FP-SBF result, but within the range of the systematic uncertainties in the various ties.

Another interesting facet of this work relates to the “fluctuation number” $`\overline{N}\overline{m}m_{\mathrm{tot}}`$, which measures the galaxy luminosity in units of the weighted mean stellar luminosity. $`\overline{N}`$ correlates tightly with stellar velocity dispersion; it also correlates with galaxy color and is independent of Galactic extinction. Interestingly, SBF distances calibrated using the properties of $`\overline{N}`$, such as those shown in Fig. 1b, can be viewed as a hybrid of SBF and FP distances, and may be more accurate than those calibrated from galaxy color alone. We plan to investigate these issues in more detail.

## 2. Cepheid Distances to Early-type Galaxies?

Cepheids occur only in spirals and other late-type, star-forming galaxies. However, the most massive virialized structures in the nearby universe (e.g., the Virgo, Fornax, and Centaurus clusters), are overwhelmingly dominated by early-type galaxies. Although some spirals appear in projection against the Virgo core, the various secondary indicators tied to the Cepheid scale indicate that these galaxies are not at the same distance as the core ellipticals (e.g., Tonry et al. 2000; Ferrarese et al. 2000; Kelson et al. 2000; Blakeslee et al. 2002). Alternatively, it may be that the secondary indicators are yielding systematically different results for the calibrating spirals and the target ellipticals.

We have an ongoing Cycle 10 WFPC2 program to calibrate the early-type galaxy distance scale via Cepheid distances to late-type galaxies that are physically associated with ellipticals. The target galaxies are the NGC 4647/NGC 4649 pair (Fig. 2a) and NGC 5128 (Cen A), an elliptical with a central dust lane and associated star formation, apparently resulting from the incursion of a gas-rich dwarf. The Cen A Cepheid observations have yielded more than 60 superb Cepheids, making this one of largest high-quality HST Cepheid data sets. At present, we are still finalizing the analysis, but Fig. 2b shows some example light curves; differential extinction within Cen A is a major issue for this program.

## 3. SNe Ia versus SBF

SBF and Type Ia supernovae (SNe Ia) studies have in the past disagreed on $`H_0`$ at the $``$20% level, which is surprising for two methods that routinely achieve 5–10% internal accuracy. The excellent resolution of HST provides an enormous advantage over ground-based data for SBF studies, and we have recently used WFPC2 to measure high-quality HST/WFPC2 SBF distances to seven early-type galaxies that have hosted well-observed SNe Ia (Ajhar et al. 2001).

The results showed excellent agreement in the relative distances, but an offset of $``$0.25 mag in zero points, which we traced to the different, and indeed dissonant, compilations of Cepheid distances used in the past for the respective zero-point calibrations of the two methods. When calibrated consistently, SBF and SNe Ia also agree in an absolute sense (Fig. 3a) and give $`H_073`$. This is the first time the agreement has been demonstrated through a direct comparison of statistically significant samples of SBF and SNe Ia galaxy distances.

## 4. The Zero Point Problem

The most pressing problem in the measurement of extragalactic distances appears to be systematic uncertainties in the zero points. We have seen that there is significant uncertainty in the zero-point tie of the early-type galaxy distance scale to Cepheids, but perhaps even greater is the uncertainty in the Cepheid zero point itself, in part due to the poorly-constrained LMC distance.

Stellar population models can be used to predict SBF magnitudes and colors for a large range of metallicities and ages (e.g., Blakeslee, Vazdekis, & Ajhar 2001b). These models reproduce the observed SBF colors and behaviors very well, but predict an SBF zero point fainter than the Cepheid-calibrated one by $`0.2\pm \mathrm{\hspace{0.17em}0.1}`$mag in $`I`$ (the only band in which SBF is directly tied to the Cepheids via spiral bulges). However, the model and empirical zero points would come into close agreement if the Cepheid scale were revised to agree with the dynamical distance to the NGC 4258 water maser (Herrnstein et al. 1999), for example, by changing the assumed LMC distance modulus from 18.5 to 18.3 mag. Further refinements of the models should provide more stringent tests of the distance scale and guide future SBF programs (Fig. 3b).

### Acknowledgments.

The projects reviewed in this paper have all been group efforts, and I thank my many distance-scale collaborators.

## References

Ajhar, E.A., Tonry, J.L., Blakeslee, J.P., Riess, A.G., Schmidt, B.P. 2001, ApJ, 559, 584

Blakeslee, J.P., Lucey, J.R., Barris, B.J., Tonry, J.L., Hudson, M.J. 2001a, MNRAS, 327, 1004

Blakeslee, J.P., et al. 2002, MNRAS, in press (astro-ph/0111183)

Blakeslee, J.P., Vazdekis, A., Ajhar, E.A. 2001b, MNRAS, 320, 193

Ferrarese, L., et al. 2000, ApJ, 529, 745

Freedman, W.L., et al. 2001, ApJ, 553, 47

Herrnstein, J. R., et al. 1999, Nature, 400, 539

Hudson, M.J., Lucey, J.R., Smith, R.J., Schlegel, D.J., Davies, R.L. 2001, MNRAS, 327, 265

Kelson, D.D., et al. 2000, ApJ, 529, 768

Tonry, J.L., Blakeslee, J.P., Ajhar, E.A., Dressler, A. 2000, ApJ, 530, 625

Tonry, J.L., et al. 2001, ApJ, 546, 681


text:
# Spectral properties of the Narrow-Line Seyfert 1 Quasar PG1211+143

## 1 Introduction

The Narrow Line quasar PG1211+143 is a very good candidate to study the broad band emission and variability properties due to a rich sample of observational data available for this object. The source exhibits particularly strong hard optical spectrum and profound steep soft X–ray emission (Elvis, Wilkes, & Tananbaum 1985). It has been frequently argued that the optical/UV and soft X–ray emission form a single Big Blue Bump component which extends across the unobserved XUV band (Bechtold et al. 1987), and strongly dominates the bolometric luminosity. However, the nature of the soft X–ray emission below $`1`$ keV is still not determined unambiguously, because any interpretation relies on extrapolating the spectrum over a decade in frequency through the unobserved XUV band.

Detailed analysis of the broad band spectra of Seyfert 1 galaxy NGC 5548 (Magdziarz et al. 1998) led to the conclusion that the optical/UV spectrum is well modeled by an accretion disk, and the hard X–rays are reproduced by standard thermal Comptonization. However, an additional component is needed to model the soft X–ray source: the most viable candidate here is another, optically thick Comptonizing medium possibly associated with the transition between the disk and the hot plasma.

Here we investigate if a similar model also applies to PG1211+143, and if this model is unique. We have analyzed the data form the ROSAT setellite, which span the period 1991-1993, ASCA observation form June 1993 and RXTE observation form August 1997. The ROSAT and ASCA spectral data were reduced using the standard software package, where for ASCA, we used data from all four detectors (SIS0/1 and GIS2/3) fitted simultaneously. For RXTE observation we use the spectrum derived from the summed Proportional Counter Array (PCA) data analyzed in a standard manner (including the background subtraction), using the ftool script rex. All analysis of the spectral X–ray data was done using the XSPEC software package, version 10.0 (Arnaud 1996). The details are given in Janiuk et al. 2001.

## 2 Results of the spectral analysis

The spectral fitting best model results are given in Tables 1 and 2. For the ASCA/ROSAT we supplement the comptt component with an underlying hard X–ray power law component characterized by its photon index, and assume fixed Galactic absorption with $`N_H=2.8\times 10^{20}`$ cm-2. In case of RXTE data we adopted the composite model consisting of a primary soft component (described as a Comptonized black body), power law continuum, and the Compton – reflected spectrum with the iron $`K\alpha `$ line included (Życki, Done & Smith 1997; 1998).

The reprocessed spectrum is parameterized by the ionization parameter $`\xi =L_X/n_er^2`$ and the reflection amplitude $`R=\mathrm{\Omega }/2\pi `$. The outer disk radius, as before, was fixed at $`10^4R_g`$ and the power law index, describing the incident flux as a function of radius was fixed at the value -3.0. However, the inclination angle and the inner radius of the reflecting disk surface were allowed to vary. The iron line energy is not a free parameter, but it depends on the ionization state of the reflecting medium. Its profile is connected with the assumed geometrical parameters: inner radius of accretion disc and viewing angle.

The amplitude of reflection is consistent with 1.0, and the reflecting area extends down to the marginally stable orbit which means that the disk is probably not disrupted. However, its surface (at least in the inner parts dominating the reflection) is considerably ionized ($`\xi 500`$), and the neutral reflection is excluded at a $`90\%`$ confidence level. For such a case, the resulting iron line energy corresponds to the domination of Fe XXII ions and its rest frame energy is $`E6.7`$ keV, while the line equivalent width is $`EW150`$ eV. The inclination angle of the disk is low, i.e. we see the disk essentially face on.

## 3 The broad band spectrum

For the purpose of a broad band, Optical/UV/X-ray model, we add a new component, which essentially has a form of standard accretion disk of Shakura & Sunyaev (1973) around a non-rotating black hole. We assume that the disk, whether it is or it is not covered by the ionized skin, radiates locally as a black body.

We assume that the emission of the standard cold disk is only detectable from radii larger than certain $`R_{ion}`$. The gravitational energy dissipated above $`R_{ion}`$ is directly re-emitted by the disk, while the energy dissipated below $`R_{ion}`$ has to provide the energy source for hard X–ray emission and Comptonizing medium. The disk parameters are therefore: the black hole mass $`M`$, the mass accretion rate $`\dot{M}`$, and $`R_{ion}`$. We assume a face-on view of the accretion disk, and the distance to the source calculated assuming $`H_o=75`$ km s-1Mpc-1.

In the upper panel of the Figure

spec2 we plot the broad band spectrum of PG1211+143, constructed for the case of the maximum value of temperature of the seed photons undergoing Comptonization. In this case there is a clear gap between the UV and soft X–ray band, which would be seen as a local minimum in the continuum around $`2\times 10^2`$ keV. The bolometric luminosity measured between 10 000 Å and 100 keV (assuming an extrapolation of the measured hard X–ray power law) is $`1.5\times 10^{45}`$ erg s-1. The mass of the black hole is $`2.2\times 10^8M_{}`$, the accretion rate is 0.4$`M_{}`$ yr-1, and $`R_{ion}=28R_g`$. $`L/L_{Edd}`$ ratio in this solution is 0.05. About 50 % of the energy is dissipated in the black body component of the disk.

In the lower panel of Figure

spec2 we show the broad band spectrum determined under the assumption that cold disk itself provides the seed photons for Comptonization, i.e. the soft photon temperature is given by the effective temperature of the inner disk. Such a spectrum forms a continuous Big Blue Bump extending from the optical/UV to soft X–ray band. The bolometric luminosity of the source is in this case higher, $`1.9\times 10^{45}`$ erg s-1. The mass of the black hole is $`1.5\times 10^8`$ $`M_{}`$, the accretion rate is 0.54 $`M_{}`$ yr-1, and $`R_{ion}=42R_g`$. $`L/L_{Edd}`$ ratio in this solution is equal to 0.10. About 30 % of the energy is directly emitted by a black body outer disk, while 60 % is in the Comptonized component and $`10`$ % in the hard X–ray power law.

Both fits presented above reproduce the broad band spectral data well. We note that all intermediate solutions are also possible and at a level of pure data analysis, we cannot reject any of the above possibilities.

## 4 Conclusions

Our spectral analysis requires the existence of a number of various plasma zones: (1) very hot plasma, responsible for the hard X–ray emission, (2) moderately hot, moderately dense plasma responsible for Comptonization, (3) moderately hot plasma responsible for reflection and (4) a cold, dense disk. It is natural to pose the question whether the components (2) and (3) are actually the same medium but parameterized differently in the spectral fits. Such interpretation clearly simplifies the geometry of the flow. The accreting stream consists simply of the outer black body disk and an inner disk which develops very optically thick ($`\tau 10`$), warm skin ($`T10^6`$K). A small fraction of energy is released in the form of hard X–ray emission above this skin. A possible schematic view of the flow is shown in Figure

fig:geom.

REFERENCES

Bechtold J., Czerny B., Elvis M., Fabbiano G., Green R.F., 1987, ApJ, 314, 699

Elvis M. Wilkes B., McDowell J.C., Green, R. F., Bechtold J., Willner S. P., Oey M. S., Polomski E., Cutri R., 1994, ApJS, 95, 1.

Janiuk A., Czerny B., Madejski G.M., 2001, ApJ, 557, 408

Magdziarz P., Blaes O.M., Zdziarski A.A., Johnson W.N., Smith D.A., 1998, MNRAS, 301, 179


text:
# Construction of raising and lowering operators for 𝒰_𝑞(𝔰𝔩(2,ℝ))

## 1 Introduction

Despite of an intensive development of the theory of quantum groups during the last decade the results concerning real noncompact forms are rather rare and far away of being complete. Even in the lowest rank cases one encounters serious difficulties (see, for example, references ). Naturally, most attention has been paid to groups with a physical interpretation . Here we concentrate on the quantum real group $`𝒰_q(𝔰𝔩(2,))`$. Even this simple example lacks, to the author’s knowledge, a consistent theory on the $`q`$-deformed level.

Let us first recall some basic facts concerning the irreducible unitary representations of the classical Lie group $`SL(2,)`$ . Let $`H,E,F`$ be the standard basis elements of $`𝔰𝔩(2,)`$. All of them are skew self-adjoint with respect to the corresponding $``$-involution in $`𝒰(𝔰𝔩(2,))`$. Set

$$J=EF,W^\pm =H\pm \mathrm{i}(E+F).$$

For any irreducible representation $`\varrho `$ of $`SL(2,)`$ in a Hilbert space $``$ the spectrum of $`\mathrm{i}\varrho (J)`$ is simple and contained in $``$. If $`\{v_n\}_{nM}`$, $`M`$, is a basis of $``$ formed by suitably normalised eigenvectors of $`\mathrm{i}\varrho (J)`$ then

$$\mathrm{i}\varrho (J)v_n=nv_n,\varrho (W^\pm )v_n=(s+1\pm n)v_{n\pm 2},$$

where $`s`$ is a fixed parameter.

Our main goal here is to propose $`q`$-deformed analogues for the generator of the compact subalgebra $`𝔰𝔬(2)`$, $`J`$, and for the raising and lowering operators $`W^\pm `$. This short contribution does not aim at all, however, to develop a complete representation theory for $`𝒰_q(𝔰𝔩(2,))`$; this task is deferred to another publication.

## 2 $`𝒰_q(𝔰𝔩(2,))`$ at root of unity

We shall consider the case when $`q`$ is a root of unity. We assume that

$$q=\mathrm{exp}(2\pi \mathrm{i}P/Q)$$

where $`Q`$ is odd, $`P\{1,2,\mathrm{},Q1\}`$ and $`P`$ and $`Q`$ are relatively prime integers. Let $`K,K^1,E,F`$ be the standard generators of the complex Hopf algebra $`𝒰_q(𝔰𝔩(2,))`$. The $``$-involutions defining the real form is prescribed by the relations

$$K^{}=K,\text{ }E^{}=q^1E,\text{ }F^{}=qF.$$

Here we prefer another set of generators generating a “slightly smaller” Hopf algebra $`𝒰`$, namely

$$X=\mathrm{i}q^1EK^1,\text{ }Y=\mathrm{i}qFK^1,\text{ }Z=K^2.$$

Thus the defining relations for $`𝒰`$ are

$`ZZ^1=Z^1Z=1,ZX=q^2XZ,ZY=q^2YZ,`$ (1)
$`q^1XYqYX={\displaystyle \frac{1}{qq^1}}(Z^21),`$
the comultiplication is defined by

$$\mathrm{\Delta }Z=ZZ,\mathrm{\Delta }X=1X+XZ,\mathrm{\Delta }Y=1Y+YZ,$$

the antipode is defined by

$$S(Z)=Z^1,S(X)=XZ^1,S(Y)=YZ^1,$$

and the counit is defined by

$$\epsilon (Z)=1,\epsilon (X)=\epsilon (Y)=0.$$

The $``$-involution defining the real form now reads

$$Z^{}=Z,X^{}=X,Y^{}=Y.$$
(2)
The irreducible representations are divided into two families (see the original paper  or the monographes ). The representations from the first family are labelled by a discrete index $`r=0,1,\mathrm{},Q1`$, the representation space is $`_r=\mathrm{span}\{v_0,v_1,\mathrm{},v_r\}`$, and the representation is defined by

$$Zv_j=\pm q^{r2j}v_j,Xv_j=q^{r2j1}[rj]_qv_{j+1},Yv_j=[j]_qv_{j1}.$$
(3)
The representations from the second family are labelled by 3 complex parameters $`\lambda ,a,b`$, the representation space is always $`=\{v_0,v_1,\mathrm{},v_{Q1}\}`$, the representation is defined by

$`Zv_j=\lambda q^{2j}v_j,`$ (4)
$`Xv_j=\mathrm{i}q^{j1}\left(ab[j]_q{\displaystyle \frac{\lambda q^{j1}\lambda ^1q^{j+1}}{qq^1}}\right)v_{j1}\text{for }j0,`$
$`Xv_0=\mathrm{i}q^1av_{Q1},`$
$`Yv_j=\mathrm{i}\lambda q^{j+1}v_{j+1}\text{for }jQ1,`$
$`Yv_{Q1}=\mathrm{i}\lambda bv_0.`$
The two families intersect for the values of parameters $`r=Q1`$ in (3) and $`a=b=0`$, $`\lambda =\pm q^{Q+1}`$ in (4).

Unfortunately, as pointed out in reference , the only representations from this list which are unitarizable with respect to the $``$-involution (2) are one-dimensional. This absolutely unsatisfactory situation of course appeals for a modified definition of $`𝒰_q(𝔰𝔩(2,))`$. A possible way out has been proposed in reference . The set of generators is augmented by introducing an additional generator $`T`$ which satisfies

$$T^2=1,\text{ }\mathrm{\Delta }T=TT,\text{ }S(T)=T,\text{ }\epsilon (T)=1.$$

The augmented Hopf algebra will be called $`𝒰^{\mathrm{}}`$. The modified $``$-involution takes the form

$$X^{}=TXT,\text{ }Y^{}=TYT,\text{ }Z^{}=TZT,\text{ }T^{}=T.$$
(5)
As shown in papers , all representations from the first family (3) are unitarizable in this modified sense and may be regarded as $`q`$-deformations of representations of $`SL(2,)`$ belonging to the discrete series.

## 3 Raising and lowering operators

Set

$$J=(qXq^1Y)Z^1=Z^1(q^1XqY).$$

One easily finds that

$$\mathrm{\Delta }J=Z^1J+J1,S(J)=ZJ,\epsilon (J)=1.$$
(6)
For the extended algebra $`𝒰^{\mathrm{}}`$ it seems to be reasonable to impose the condition

$$JT=TJ,$$

i.e. $`J^{}=J`$.

Note that conversely one can express

$$X=\frac{1}{q^2q^2}(qJZq^1ZJ),Y=\frac{1}{q^2q^2}(q^1JZqZJ).$$

It is also straightforward to derive the relations between $`J`$ and $`Z`$. The second and third relations in (1) give

$$Z^2J(q^2+q^2)ZJZ+JZ^2=0$$
(7)
while the last one implies

$$(q^2+1+q^2)ZJ^2ZJZJZJZ^2JZJZJ=(q+q^1)^2(Z^21).$$
(8)
All the above manipulations are reversible and this shows

###### Proposition 1.

The Hopf algebra $`𝒰`$ is isomorphic to the Hopf algebra $`𝒰^{}`$ generated by $`Z`$ and $`J`$ with the defining relations (7), (8), (6) and the known relations for $`Z`$ ($`ZZ^1=Z^1Z=1`$, $`\mathrm{\Delta }Z=ZZ`$, $`S(Z)=Z^1`$, $`\epsilon (Z)=1`$).

The quantum numbers $`[x]_q=(q^xq^x)/(qq^1)`$ are supposed to be defined for any $`x`$ for we set, by definition, $`q^x=\mathrm{exp}(2\pi \mathrm{i}(P/Q)x)`$. The main result of this contribution is the following identity.

###### Proposition 2.

For all $`x`$ it holds true that

$`Z\left(J[x+2]_q\right)\left(J[x]_q\right)\left(J[x2]_q\right)Z`$ (9)
$`=\left(\left(J[x]_q\right)Z\left(J[x]_q\right)Z[2]_{q}^{}{}_{}{}^{2}\right)\left(J[x]_q\right),`$
###### Proof.

To prove (9) we expend the both sides in the variable $`y=q^x`$,

$$\mathrm{LHS}\mathrm{RHS}=c_2y^2+c_1y^1+c_0+c_1y+c_2y^2,$$

and we find that

$`(qq^1)^2c_2=(qq^1)^2c_2=JZ^2+(q^2+q^2)ZJZZ^2J,`$
$`(qq^1)c_1=(qq^1)c_1`$
$`=(q^2+1+q^2)ZJ^2ZJZJZJZ^2JZJZJ[2]_{q}^{}{}_{}{}^{2}(Z^21),`$
and

$`(qq^1)^2c_0=`$ $`\mathrm{\hspace{0.17em}2}(JZ^2(q^2+q^2)ZJZ+Z^2J)`$
$`+(qq^1)^2\left(ZJ^3Z[2]_{q}^{}{}_{}{}^{2}ZJZJZJZJ+[2]_{q}^{}{}_{}{}^{2}J\right).`$
Thus the coefficients $`c_{\pm 2}`$, $`c_{\pm 1}`$ vanish owing to (7) and (8), the coefficient $`c_0`$ vanishes if it is true in the particular case $`x=0`$. The following lemma concludes the proof. ∎

###### Lemma 3.

It holds true that

$$Z\left(J^2(q+q^1)^2\right)JZ=\left(JZJZ(q+q^1)^2\right)J.$$

###### Proof.

To prove the particular case of (9) with $`x=0`$ set

$$V=ZJ^3Z[2]_{q}^{}{}_{}{}^{2}ZJZ+[2]_{q}^{}{}_{}{}^{2}JJZJZJ.$$

Then it holds

$`ZV+VZ=`$
$`\left((q^2+1+q^2)ZJ^2ZJZ^2JZJZJJZJZ+(q+q^1)^2(1Z^2)\right)`$
$`\times JZ+ZJ`$
$`\times \left((q^2+1+q^2)ZJ^2JZ^2JZJZJJZJZ+(q+q^1)^2(1Z^2)\right)`$
$`+(JZ^2(q^2+q^2)ZJZ+Z^2J)J^2Z`$
$`+ZJ^2(JZ^2(q^2+q^2)ZJZ+Z^2J)`$
and thus $`ZV+VZ=0`$. This clearly implies

$$Z^kV=(1)^kVZ,\text{for all }k.$$

In particular, $`Z^QV=VZ^Q`$. But $`Z^Q`$ belongs to the centre of the algebra $`𝒰`$ and therefore $`Z^QV=VZ^Q`$. Hence $`Z^QV=0`$ and consequently $`V=0`$. ∎

Some straightforward consequences follow from identity (9). Consider an irreducible representation of $`𝒰`$ whose dimension equals $`d`$. If $`v`$ is an eigenvector of $`J`$ with an eigenvalue $`[x]_q`$ then the vectors

$$\left(J[x]_q\right)\left(J[x2]_q\right)Zv,\left(J[x]_q\right)\left(J[x+2]_q\right)Zv,$$
(10)
either vanish or are eigenvectors of $`J`$ corresponding to the eigenvalues $`[x+2]_q`$ and $`[x2]_q`$, respectively. The spectrum of $`J`$ is of the form

$$[x]_q,[x+2]_q,\mathrm{},[x+2d2]_q,$$

and hence the operator $`J`$ has similar properties (up to the factor $`\mathrm{i}`$) as its nondeformed counterpart. For example, in the case of the representations from the first family (3) the operator $`J`$ is known to have the spectrum (10) with $`x=Qd+1`$ . Furthermore, the operators $`(J[x]_q)(J[x2]_q)Z`$ and $`(J[x]_q)(J[x+2]_q)Z`$ play the role of raising and lowering operators, respectively. Finally, the matrix $`Z`$, if expressed in the corresponding eigenbasis of $`J`$, is tridiagonal in the case of the first family of representations (3) or tridiagonal with the cyclic convention modulo $`Q`$ in the case of the second family (4).

## Acknowledgements

The research was partially supported by Grant GACR 201/01/01308.


text:
# Unified Scaling Law for Earthquakes

\[

## Abstract

We show that the distribution of waiting times between earthquakes occurring in California obeys a simple unified scaling law valid from tens of seconds to tens of years, see Eq. (1) and Fig. 4. The short time clustering, commonly referred to as aftershocks, is nothing but the short time limit of the general hierarchical properties of earthquakes. There is no unique operational way of distinguishing between main shocks and aftershocks. In the unified law, the Gutenberg-Richter $`b`$-value, the exponent $`1`$ of the Omori law for aftershocks, and the fractal dimension $`d_f`$ of earthquakes appear as critical indices.

\]

Earthquakes are a complicated spatio-temporal phenomenon. The number of earthquakes with a magnitude $`M>m`$ is given by the Gutenberg-Richter law . In addition to the regularity in the rate of occurrence, earthquakes display a complex spatio-temporal behavior . The spatial distribution of epicentres is fractal and they occur on a fractal-like structure of faults . Short-range temporal correlations between earthquakes are expressed by Omori’s law , which states that immediately following a main earthquake there is a sequence of aftershocks whose frequency decays with time as $`T^\alpha ,\alpha 1`$. This has led to the commonly held belief that aftershocks are caused by a different relaxation mechanism than the main shocks.

The observed temporal complex behavior is obviously of dynamical origin. However, the statistics of earthquakes as well as the geometrical fractal structure displayed by the faults and by the spatial distribution of epicenters is also a result of a dynamical process and one might speculate whether it is possible to unify these observations.

We propose a unified scaling law for the waiting times between earthquakes, expressing a hierarchical organization in time, space, and magnitude. There is a correlated regime where the distribution of waiting times between earthquakes is a power-law, $`T^\alpha ,\alpha 1`$ and an uncorrelated regime. However, the waiting time interval for the crossover between the two regimes for earthquakes larger than a given magnitude depends on the area and magnitude under consideration.

An earthquake catalogue covering the period $`19842000`$ in a region of California spanning 20N - 45N latitude and 100W - 125W longitude was analyzed . The total number of recorded earthquakes in the catalogue is $`335,076`$. The number of earthquakes $`N(M>m)`$ with magnitude larger than $`m`$ is given by the Gutenberg-Richter law  $`\mathrm{log}_{10}N(M>m)bm`$, $`b0.95`$ see Fig. 1.

The spatio-temporal analysis was carried out as follows. We covered the region with a grid with cells of size $`L\times L`$, see Fig. 2 and defined the waiting time $`T`$ as the time interval between the beginning of two successive earthquakes. We then measured $`P_{S,L}(T)`$, the distribution of waiting times $`T`$, between earthquakes occurring within range L whose magnitudes are greater than $`m=\mathrm{log}(S)`$ . \[

\] Figure 3 shows the resulting set of curves $`P_{S,L}(T)`$, for time scales ranging from seconds to $`16`$ years, for several values of $`S`$ and $`L`$, plotted on double logarithmic scale. Obviously, the curves differ widely. Some general trends can be seen, however. There is a linear regime, indicating a power-law distribution, extending up to a cut-off indicating an upper limit of the waiting time. For fixed cell size $`L`$ and increasing cut-off $`S`$ (or $`m`$), the range of the power-law regime increases. For fixed cut-off $`S`$ and increasing cell size $`L`$, the range of the power-law regime decreases .

In Fig. 4, the curves are re-plotted in terms of re-scaled coordinates. The $`x`$-axis is chosen as $`x=TS^bL^{d_f}`$, and the $`y`$-axis represents $`y=T^\alpha P_{S,L}(T)`$. The rescaling causes a shift of the curves in Fig. 3 that depends on $`L`$ and $`S`$. For suitable choice of the interval exponent $`\alpha `$, the magnitude exponent $`b`$, and the spatial dimension $`d_f`$, all the data collapse nicely onto a single well-defined curve $`f(x)`$, that is,

$$T^\alpha P_{S,L}(T)=f(TS^bL^{d_f}).$$
(1)
This equation expresses the unified scaling law for earthquakes. The function $`f(x)`$ consists of a constant part and a decaying part, separated by a sharp kink. The constant part corresponds to the linear, power-law part in Fig. 3 since we have multiplied $`P_{S,L}(T)`$ with $`T^\alpha `$. Any deviation from power-law behavior would show up dramatically in this type of plot. Nevertheless, the function is approximately constant over $`8`$ orders of magnitude! The rapidly decaying part is consistent with an exponential decaying function implying an uncorrelated regime for large values of $`x`$. This is indeed what one would expect on physical ground: earthquakes that are separated by large enough distances or long waiting times will be uncorrelated.

The index $`\alpha 1`$ can be identified as the Omori-law exponent for aftershocks, $`b1`$ is the $`b`$-value in the Gutenberg-Richter law, and $`d_f1.2`$ describes the $`2d`$ fractal dimension of the location of epi-centers projected onto the surface of the Earth.

The data collapse implies that the waiting time distribution depends on $`T,S`$, and $`L`$ only through the variable $`x`$. Only critical processes exhibit this type of data collapse, known as scaling in critical phenomena , so our analysis demonstrates that earthquakes are a self-organized critical (SOC) phenomenon , as had been anticipated from the existence of the Gutenberg-Richter law . The data collapse shows that there is no separate relaxation mechanism for aftershocks. The three exponents $`\alpha ,b`$, and $`d_f`$ characterizing earthquakes emerge as critical indices in the unified law. The estimated uncertainty in the values of the critical indices is less than 0.2. Whether the critical exponents vary with region and maybe even with time is an interesting question that is outside the scope of the present letter but we urge further studies in that direction.

Depending on the value of scaling argument $`x`$, and thus the chosen values of $`L`$ and $`m`$ (or $`S`$), two successive earthquakes will either be correlated, for $`x`$ small, (i.e. to the left of the kink in Fig. 4), or uncorrelated, for $`x`$ large (i.e. to the right of the kink in Fig. 4).

Depending on the length scale $`L`$ of observations, and the magnitude magnitude $`m`$ (or $`S`$) chosen, the correlated Omori $`T^\alpha `$ regime may range from seconds to tens of years (and probably much longer if data were available). If the earthquakes are correlated they may be interpreted as belonging to an aftershock sequence. If they are uncorrelated, they may be interpreted as main events. However, this depends on $`L`$ and $`m`$ through the variable $`x=TS^bL^{d_f}`$ and has no absolute meaning. Therefore, there is no unique way of characterizing earthquakes as aftershocks or main events.

To summarize, the short time correlations given by Omori’s law is just the short time limit of a general hierarchical scaling phenomenon occurring at all accessible time scales. Amazingly, the statistics of aftershocks occurring within minutes of an earthquake can be simply related to the statistics of earthquakes separated by tens of years!

One may think of the value of $`m`$ (or $`S`$) at the kink as being a ’characteristic’ magnitude of earthquakes. However, since the position of the kink is a function of $`x`$, not $`m`$ (or $`S`$), one cannot specify the magnitude of the ’characteristic’ earthquake without at the same time specifying a time scale and an area. However, there is no special time scale that can play any absolute role for the dynamics of earthquakes, limited at the upper end by the time scale of tectonic plate motion, and at the lower end by the duration of earthquakes.

How should one physically understand the fundamental law, Eq. (1)? Let us first discuss the meaning of the scaling variable $`x=TS^bL^{d_f}`$. The quantity $`S^bL^{d_f}`$ appearing in it’s definition is a measure of the average number of earthquakes per time unit with magnitude greater than $`m=log_{10}(S)`$ occurring within the range $`L`$. Thus, $`x`$ is a measure of the average number of such earthquakes occurring within a time interval $`T`$. The law states that the distribution of waiting times depends only on this number. When this number exceeds a well-defined value (the position of the kink in Fig. 4), the earthquakes sharply become less correlated.

Think of earthquakes being generated by ’processes’, each producing a sequence of correlated earthquakes with a $`T^\alpha `$ distribution. These processes correspond to a sequence of avalanches in self-organized critical models of complex phenomena. Visually, one might think of the processes as the activity associated with dynamically changing fault segment patterns. The law indicates that the crust operates in the true SOC slow-driving regime  where the individual processes (avalanches) do not overlap. Because of the non-zero driving rate, several spatially separated processes are active simultaneously. The kink on the $`f`$-curve indicates the point where one crosses into the regime where spatially independent earthquakes, belonging to different processes, are sampled within a window spanned by $`T`$ and $`L`$. For small enough $`L`$ and $`T`$ only a single, correlated process is sampled.

Acknowledgements. We would like to thank C. Scholz for advising us on earthquake catalogues. P.B. acknowledge helpful discussions with M. Paczuski. Correspondence and requests for materials should be addressed to P.B. (email: bak@alf.nbi.dk) or K.C. (e-mail: k.christensen@ic.ac.uk).


text:
# Principal Component Analysis for Spectral Indices of Stellar Populations

## 1 Introduction

To understand how galaxies are formed and evolve we need to investigate their stellar populations and to derive their main parameters, such as metallicity and age. This study plays an important role in our understanding of galaxy properties (Kennicutt 1998; Maraston 1998; Padoan et al. 1997). During the last decade particularly using population synthesis method has performed determination of galaxy properties. It has been used extensively by many authors for all kinds of galaxies (Leitherer et al. 1996). However, key parameters in interpretation of the observed properties are the metallicity and the age. The problem is the degeneracy of the effects from variations in age and in metallicity, even in the simplest unresolved stellar systems their effects are very difficult to separate (Vazdekis et al. 1997). It makes the determination of age and metallicity uncertain. To determine the age and metallicity more accurately, strong efforts have been devoted to select some spectral features that are more sensitive to age, and others are more sensitive to metallicity. In 1994, G. Worthey developed a very efficient method to investigate the age and metallicity effects in the integrated light of stellar populations. In his method, suppose an index varies by an amount $`\mathrm{\Delta }I`$, it can be explain $`\mathrm{\Delta }I`$ by either a pure metallicity variation or a pure age variation. The Z sensitivity parameter is the ratio of the percentage change in age to the percentage change in Z, with larger numbers indicating greater metallicity sensitivity (Worthey 1994). Using this method, some metallicity sensitive spectral indices are found, but only two age sensitive spectral indices are found.

In this paper, we apply a different statistical technique, principal component analysis (PCA) to simple stellar population sample, to investigate on the reliability of the relations between spectral indices and age. We want to extract other spectral indices that can be used to determine the age of stellar population. The samples of simple stellar population (SSP) spectra come from Bressan et al. (1996) and Bruzual & Charlot (1996). For each spectra, 21 spectral indices in the Lick/IDS system be measured first (Trager et al. 1998), then we use PCA method to find some spectral indices that are more sensitive to the age. We leave the synthetic galaxy spectral indices with different star formation history and observational galaxy spectral indices to a future paper. However we already notice an encouraging good resemblance between the results of PCA and those one gets from observation (Worthey 1994).

The organization of the rest of this paper is as follows. In section 2 we present PCA method. In section 3 we describe the definitions of spectral indices, the synthetic models we have used. In section 4 we analyze a sample of SSP spectra with different ages, find the first principal components (PCs), discuss the significance of different spectral indices to determine the age, and how different age sampling affects the results. In section 5, we illustrate that the uncertainties will not affect our results, and compare our results with previously published values. The conclusions are summarized in section 6.

## 2 Principal Component Analysis

Principal component analysis is known to be powerful technique for unveiling correlation between variables in a data set and for determining the intrinsic dimensionality of a parameter space (Connolly & Szalay 1999). It is an orthogonal transformation that allows the building of more compact linear combinations of data that are optimal with respect to the mean square error criterion. The algorithm calculates a new base with the minimum set of orthogonal axes that describes the observed variance of the data. Finding the new axes is an eigenvalue formation matrix between the old base and the new base. The new orthogonal basis is composed of vectors called principal components (PCs). The corresponding eigenvalues represent the variances of the parameters in the new bases (Francis & Wills 1999).

The formulation of standard PCA is as follows (see e.g. Murtagh & Heck 1987). Consider a set of $`n`$ objects ($`i=1,n`$), each with $`m`$ parameters ($`j=1,m`$). If $`x_{ij}`$ are the original measurements, the starting point in PCA method is an observation matrix $`X=(x_{ij})_{n\times m}`$ in which row vectors list observations characterizing an object.

We then construct a covariance matrix $`C=(c_{ij})_{m\times m}`$, where

$$c_{jk}=\frac{1}{n1}\underset{i=1}{\overset{n}{}}(x_{ij}\overline{x}_j)(x_{ik}\overline{x}_k)1j,km$$
(1)
where $`\overline{x}_j`$ and $`\overline{x}_k`$ are the mean of column vector.

It can be shown that the axis along which the variance is maximal is the eigenvector $`𝒆_\mathrm{𝟏}`$ of the matrix equation

$$C𝒆_\mathrm{𝟏}=l_1𝒆_\mathrm{𝟏}\mathbf{,}$$
(2)
where $`l_1`$ is the largest eigenvalue of matrix C, which is in fact the variance along the new axis. The other principal axes and eigenvectors obey similar equations. It is convenient to sort them in decreasing order, and to quantify the fractional variance by $`l_j/_jl_j`$. The matrix of all the eigenvectors forms a new set of orthogonal axes that are ideally suitable to a description of the data set. If the main variance of the data set lies in a small dimensional space, then one can get a good visualization of it by plotting the data projected on the first few eigenvectors (PCs). The projection of a vector $`𝒙`$ on the eigenvector $`𝒆_𝒋`$ being $`𝒙\mathbf{}𝒆_𝒋`$, where $`𝒙`$ is a row vector of $`X`$. The fractional variance of the first PCs tells us to what extent the data lie in a given low dimensional space (Ronen et al. 1999).

In what follows, we will apply this technique to different spectral index samples. To avoid confusion we denote eigenvectors by PC1, PC2 etc., and projections on the new axes by pc1, pc2 etc. For the $`k`$th object, $`pc_j=𝒙𝒆_𝒋=e_{j1}x_{k1}+\mathrm{}+e_{ji}x_{ki}+\mathrm{}+e_{jm}x_{km}`$.

## 3 Indices of Simple Stellar Population

Among the libraries of spectral indices, the most extensive one is created by Faber and coworkers at Lick Observatory, the Lick/IDS System (Faber et al. 1985; Worthey et al. 1994; Trager et al. 1998). In the Lick/IDS system, absorption-line strengths are described by “indices,” in which a central “feature” band-pass is flanked to the blue and red by “pseudo-continuum” band-passes. It is fully described in Worthey et al. (1994) and Trager et al. (1998). Table 2111This table is available electronically from the Astrophysical Data Center (http://adc.gsfc.nasa.gov/adc.html). of Trager et al. (1998) presents the band-passes of these 21 Lick/IDS absorption line indices and the features measured by these indices.

Simple stellar population is defined as a single generation of coeval stars with fixed parameters such as metallicity, initial mass function, etc (Buzzoni 1997). To find some spectral indices that can be used to determine the age of stellar population, as the first step, we apply PCA method to SSP only. To assess how model- dependent our conclusions are, we use the spectral indices of SSP form two different synthetic models.

The first SSP sample is the one used by Tantalo et al. (1996). Detailed descriptions of the SSP model can be found in Bressan et al. (1994), Silva (1995), and Tantalo et al. (1996). These SSPs extend the ZAMS to 0.15 $`M_{}`$. The rate of mass loss during the RGB and TP-AGB phases is according to the Reimers (1975) relationship with $`\eta =0.45`$. They are computed with the new library of stellar spectra described in Silva (1995), and Tantalo et al. (1996). These SSPs are shortly referred to as SIL-SSP. The spectral indices of this sample are calculated by means of the empirical calibrations of Worthey (1992) and Worthey et al. (1994).

The second one is isochronal synthesis model of Bruzual & Charlot (1996). It has extended the Bruzual & Charlot (1993) evolutionary population synthesis model. The updated version provides the evolution of the spectrophotometric properties for a wide range of stellar metallicity. The stellar evolution tracks are from Padova group (Bressan et al. 1993). The stellar spectra library is from Lejeune et al. (1997,1998). These SSPs are shortly referred to as BC-SSP.

The spectral indices in Lick/IDS system for SIL-SSP and BC-SSP are published in the AAS CD-ROM Series, Vol. 7. The technique to calculate spectral indices of SSP is amply described in Worthey (1994), Bressan et al.(1996), Tantalo et al. (1998) and Longhetti et al. (1998).

## 4 PCA to SSP

To find some spectral indices that are more sensitive to the age, we will apply PCA method to SSPs in this section.

### 4.1 PCA to SIL-SSP

Given an IMF and metallicity, evolutionary synthesis model provides SSP’s spectra as a function of age in the sample of SIL-SSP. Each data set contains about 50 SSPs, corresponding to 50 time steps from 0.004 to 20 Gyr. Each SSP has 1206 wavelength points in the range from 91 Å to 100 $`\mu m`$ (Leitherer et al. 1996). Because the spectral indices of CN1, CN2 and TiO2 are equal to 0.0 or very little for young stellar population, and there are some uncertainties at high $`Z`$ SSP, so we perform PCA on a sample of SSP with age 1–15 Gyr, and Z=0.0004 to 0.05.

For our first study, we performed PCA on the SIL-SSP with $`Z=Z_{}=0.02`$. The results are shown in Table 1. Column (2)–(3) show the first 2 principal components (PCs) out of a total of 21 principal components. The 3rd row in Table 1 gives the variances (eigenvalues) of the data along the direction of the corresponding principal component. By convention, the principal components are given in order of their contribution to the total variance, it is given as ‘Proportion’ in the 4th row. The columns of numbers for each principal component represent the weights assigned to each input variable. Thus pc1 $`=0.02\times x_1+0.01\times x_2+0.14\times x_3+0.49\times x_4+0.49\times x_5`$ …, where $`x_1`$, $`x_2`$, $`x_3`$, $`x_4`$, and $`x_5`$ …, are the values of the normalized variables corresponding to CN1, CN2, Ca4227 ,G4300 and Fe4383, etc. If the value in Column 2 of Table 1 is large, it suggests that the corresponding spectral index is important to PC1. The first principal component is elongated with variance 5.11, and accounts for about most the total variance. It is therefore likely to be highly significant. In SIL-SSP sample, only age is variable, so PC1 shows the age information, it correlates with ages. We can use PC1 to determine the age of a SSP.

This result shows that most of these spectral indices correlate with PC1, but the correlation involving G4300, Fe4383, C24668, H$`\beta `$, and Mg$`b`$ are very strong. So these spectral indices can be used to determine the age of stellar system. PC2 accounts for only 1% of the variance in this data. Some spectral indices appear to contribute to both PC1 and PC2, but the contribution of PC2 is very little, therefore any spectral indices present in PC2 represent only fine tuning to the main variation by which PC1 represents. PC2 is not a significant component.

To inspect how sensitive the results of PCA method are to the different metallicity, we performed PCA on the SIL-SSP with other metallicity. The results of performing PCA on these SSPs are shown in the last 8 columns of Table 1. There are some distinct characters in the results. First, we find that the spectral indices G4300, Fe4383, C24668, H$`\beta `$, and Mg$`b`$ are very strong correlated with PC1 also. Second, the contribution from PC1 is much more important than that from PC2. Third, we find that the spectral indices contributing to the PCs are different from the SSPs with different metallicity. Mg$`b`$ and C24668 contributing to low metallicity SSP are small, they cannot fit to determine the age for low metallicity SSP but for high metallicity SSP. On the contrary, G4300 and H$`\beta `$ are more fit to determine the age for low metallicity SSP. Actually, the sodium index Na$`D`$ is also an important index, with respect to the PC1. We did not select it due to the considerable effect of the interstellar sodium line as previously reported(Worthey 1994)

### 4.2 PCA to BC-SSP

There exist a number of SSP models that are synthesized with different approaches. It is important to check the sensitivity of our results to the SSP models adopted. The Bruzual & Charlot’s model is suitable for the comparison, because it uses a similar technique of “isochronal synthesis” to predict the spectral evolution of stellar population. Each data set (metallicity fix) contains about 221 SSPs, corresponding to 221 time steps from 0 to 20 Gyr. Each SSP has 1206 wavelength points in the range from 5 Å to 100 $`\mu m`$ (Bruzual & Charlot 1996). Similar to SIL-SSP, we performed PCA method to the stellar population (1–15 Gyr) with metallicity Z=0.0004 to 0.05. The results are shown in Table 2.

In Table 2, the first principal component also accounts for about most the total variance, more than 98%. PC2 can be ignored, it maybe caused some uncertainties. This result is consistent with the principle of PCA method: One of important application of PCA is known to be powerful technique for determining the intrinsic dimensionality of a parameter space. Given an IMF and metallicity, SSP has one parameter, age. This parameter can be expressed by PC1. From Table 1 and Table 2, we find the results from different SSP are very similar, it can be explained that our results are not sensitive to the SSP models adopted. The spectral indices, G4300, Fe4383, C24668, H$`\beta `$, and Mg$`b`$, can be used to determine the age of stellar population.

### 4.3 PCA to the SSPs with different metallicity

To extract an objective information on the age-metallicity sensitivity of the narrow band indices, we apply PCA method to the whole set of SSPs. As a result, we find a little information for metallicity. We think it can be understood easily, because there are only five kind metallicities for the sample of SSP, the metallicity point is less than the kind of spectral indices. So we can not obtain more metallicity information from this statistical technique. In the future paper, we will observe some spectra of globular cluster with different metallicity, and perform a principal component analysis method to them, It maybe help us to get the metallicity sensitive indices.

To overcome this difficulty, we have selected 3 spectral indices from Worthey (1994). There are H$`\beta `$ which is age sensitive, Fe5015 which is metallicity sensitive and Ca4227 which isn’t age or metallicity sensitive. We apply PCA method to these 3 indices, with same age and different metallicities(5 kind metallicities). The results are shown in Column (3)–(8) of Table 3 for SIL-SSPs, and in Column (9)–(14) for BC-SSPs. The age of SSP is $`T=10^9`$ yr, $`10^{10}`$ yr, and $`1.5\times 10^9`$ yr. The 1st row in Table 3 gives the age of the SSP, the meaning of next 3 rows same as Table 1. The numbers in the last 3 rows represent the weights assigned to each input variable. If the value in these columns is large, it suggests that the corresponding spectral index is important. From this table, we find the first principal component accounts for about most the total variance. This is therefore to be highly significant. In this SSP sample, since the age of SSP is same, only metallicity is variable, so PC1 must show the metallicity information. In the table, the value of the normalized variable corresponding to Fe5015 is large, it suggests that the corresponding spectral index, Fe5015, is important to metallicity. The value of H$`\beta `$ is small, it seems that H$`\beta `$ is not important to metallicity. This result is the same as the result of Worthey (1994). From this analysis, we know PCA is a very efficient method, it can be used to select some spectral indices that will help us to determine the age or metallicity of SSP, if the age and metallicity points are enough.

### 4.4 PC1 and Index with Age

For each SSP sample (metallicity fix), we can obtain $`\stackrel{}{e_j}`$ from PCA method. Using $`pc_j=\stackrel{}{x}\stackrel{}{e_j}`$, we can calculate the projection of a vector $`\stackrel{}{x}`$ on the eigenvector. We plot the projection on the PC1 (pc1) versus age for the metallicity range from Z=0.0004 to Z=0.05 in Fig. 1. Each line pattern represents a different metallicity, as indicated inside the frame. In Fig. 1(a), we plot the pc1 versus age for SIL-SSP. In Fig. 1(b), we plot for BC-SSP. In these figures, the following remarks can be made. (1) It is apparent that there is a uniform tendency for the values of pc1 to become larger as the metallicity increases from $`Z=0.0004`$ to $`Z=0.05`$. (2) These mainly are monotonically increasing curve with age, clearly. Therefore, once we know the metallicity, we can use this pc1 to determine the age of stellar population. (3) Some small non-monotonic irregularities in the curve are probably due to abrupt transition of stars to different evolutionary stages. They may be real, but they may also occur if the model does not contain enough evolutionary tracks in a certain range. In any case, we shall see that these irregularities will be smoothed out when considering longer and more realistic star formation histories instead of the instantaneous bursts considered here.

As an example, we plot the spectral indices with age for SIL-SSP with Z=0.02 in figure 2. We find that spectral indices become stronger with increasing age, while Balmer line, H$`\beta `$, become weaker. The change of G4300, Fe4383, C24668, H$`\beta `$, and Mg$`b`$ with age are very strong, especially for G4300 and H$`\beta `$. These results, from Fig. 1 and Fig. 2, consistent with the results from PCA method. So these spectral indices can be used to determine the age of stellar system. The results suggested that PCA is efficacious method for find the age sensitive indices.

## 5 Discussion

In this paper, as the first step, we only performed PCA to SSP sample. There are two reasons. Firstly, SSP is simple and reasonably well understood, so it is important to see what one can learn by using this simplest assumption, and then check whether more complex star formation histories give qualitatively similar conclusions. This is a common approach often taken in the evolutionary population synthesis models for galaxies (Vazdekis et al. 1997). Secondly, the selected spectral indices in this paper can be used to determine the age of star cluster.

The synthesis spectral indices are computed from the SSP of various ages and metallicity, under some fitting functions. The uncertainties in the fitting function, the stellar evolution input and the evolutionary population synthesis computational procedure may be affecting our result. But for the result of Maraston et al. (1999), the discrepancies introduced by different fitting function are relatively small. Form the result state above , the difference between the SIL-SSP and BC-SSP is very small, so the uncertainties due to the use of different evolutionary population synthesis procedures and stellar evolution input can be ignored.

Compared with the previous methods, PCA method provides a more objective and informative alternative to diagnostics by individual spectral lines, the PCs represent all the spectral features simultaneously. Same as Worthey (1994), we also find G4300 and H$`\beta `$ are age sensitive indices that are useful in determining age. In addition, we also find some other indices that can be used to determine the age of SSP with different metallicities.

## 6 Conclusions

In order to select spectral indices that more sensitive to age than metallicity, PCA method has been applied to the Lick indices in the population model of SIL-SSP and BC-SSP. It is the first time this method is used for selecting the age sensitive spectral indices. Below, we summarize our main conclusions.

1. Using PCA method, we find some spectral indices, such as G4300, H$`\beta `$, C24668 and Mg$`b`$, which can be used to determine the age of stellar population (1–15 Gyr). These spectral indices will help us to determine the age of stellar population more accurate.
2. Important differences are found between the age sensitive indices at different metallicities. For example, in low metallicity C24668 and Mg$`b`$ are not fit to determine the age of stellar population, but can be used to determine the age for high metallicity stellar population.
3. To understand how sensitive our method is to different stellar population synthesis models and metallicities, we have compared the results obtained with SIL-SSP and BC-SSP with 5 kind of metallicities. We find that although the precise values of PCs are difference, the general trend is very similar. We can conclude that the uncertainties in the evolutionary population synthesis models and metallicities can be ignored.

## Acknowledgments

We deeply thank the referee, Dr. A. Bressan, for so carefully reading an early draft of this paper and for making some very useful comments and constructive suggestions, which improved the manuscript. We are indebted to Dr. S. Ronen and Dr. G. Worthey to have made available their program. We are grateful to Dr. A. Bressan for providing us with a grid of spectral indices for simple stellar population. This work is supported by the Chinese National Natural Science Foundation (CNNSF) and Chinese National Pandeng Project.