text:
# Orbit bifurcations and wavefunction autocorrelations11footnote 1Short title: Bifurcations and autocorrelations

## 1 Introduction

The local statistical properties of quantum eigenfunctions $`\psi _n(𝐪)`$ in classically chaotic systems are often modelled in the semiclassical limit by treating the eigenfunctions as random superpositions of plane waves (Berry 1977). This allows, for example, for value distributions and autocorrelation functions to be computed. However, it ignores contributions from the scarring of the eigenfunctions by classical periodic orbits (Heller 1984, Kaplan 1999). Understanding these contributions has been one of the central goals in quantum chaology. Our purpose here is to study the influence of bifurcating periodic orbits on the autocorrelation of eigenfunctions.

The semiclassical theory describing the scarring of wavefunctions by isolated and unstable periodic orbits was developed by Bogomolny (1988), who derived an explicit formula for the contribution from a given such orbit to $`|\psi _n(𝐪)|^2`$, averaged locally with respect to position and energy. Berry (1989) extended this approach to the Wigner functions associated with energy eigenstates, again averaged locally with respect to the energy (position averaging is unnecessary in this case). Berry’s formula also describes the influence of isolated and unstable periodic orbits on the autocorrelation of the eigenfunctions, because the Wigner function and the autocorrelation function $`\psi _n(𝐪+𝐋)\psi _n^{}(𝐪)`$ are fourier transforms of each other (explicit formulae my be found in Li & Rouben 2001). It was shown in Keating (1991) that Bogomolny’s and Berry’s formulae are exact, rather than semiclassical approximations, for quantum cat maps.

The extension of Bogomolny’s theory to bifurcating periodic orbits was made in Keating & Prado (2001). There it was shown that orbits undergoing bifurcation give rise to scars in $`|\psi _n(𝐪)|^2`$ (averaged locally with respect to position and energy) with a width that scales as $`\mathrm{}^\omega `$ and an amplitude that scales as $`\mathrm{}^\alpha `$ when $`\mathrm{}0`$, where $`\omega `$ and $`\alpha `$ are bifurcation dependent scar exponents whose values are positive and less than or equal to those ($`\omega =\alpha =1/2`$) corresponding to isolated and unstable orbits. For most of the generic bifurcations $`\omega <1/2`$ and $`\alpha <1/2`$, and so bifurcations may be said to give rise to superscars. This mirrors the dominant influence they have previously been shown to exert on spectral fluctuation statistics (Berry et al. 1998, Berry et al. 2000).

The question we address here is how bifurcations influence the autocorrelation of the eigenfunctions. Our main result is that the length-scale of the correlations around a bifurcating orbit scale as $`\mathrm{}^{1\alpha }`$ when $`\mathrm{}0`$, where $`\alpha `$ is the scar amplitude exponent. By comparison, the random wave model predicts correlations with a length-scale of the order of the de Broglie wavelength (i.e. of the order of $`\mathrm{}`$). We illustrate the influence of scarring on the quantum autocorrelation function with numerical computations for a family of nonlinearly perturbed cat maps (Basílio de Matos & Ozorio de Almeida 1995; Boasman & Keating 1995).

The bifurcation of periodic orbits is a characteristic phenomenon in systems whose phase space is mixed. It has been conjectured by Percival (1973) that the energy eigenstates of mixed systems separate semiclassically into those associated with stable islands (regular) and those associated with the chaotic sea (irregular). The random wave model has recently been extended to irregular states in mixed systems by Bäcker & Schubert (2002) (BS1 in the references). The autocorrelation of these states was investigated in this context in Bäcker & Schubert (2002) (BS2 in the references). Our results, concerning the length scale of the correlations, are expected to apply to the autocorrelation of irregular eigenstates in the vicinity of superscars.

## 2 Autocorrelation Formulae

Our purpose here is to derive semiclassical autocorrelation formulae for bifurcating periodic orbits which generalize those obtained by Berry (Berry 1989) (see also Li & Rouben 2001) for periodic orbits far from bifurcation.

The autocorrelation function for the eigenfunction $`\psi _n(𝐪)`$ of a given quantum Hamiltonian corresponding to an energy level $`E_n`$ is

$$c_n(𝐋)=\psi _n(𝐪+𝐋)\psi _n^{}(𝐪),$$
(1)
where $`\mathrm{}`$ denotes an average over position $`𝐪`$. We shall be particularly concerned with the case when the autocorrelation function is averaged over a group of eigenstates with energy levels near to $`E`$, that is with

$$C(𝐋,E,\epsilon )=\underset{n}{}c_n(𝐋)\delta _\epsilon \left(EE_n\right),$$
(2)
where $`\delta _\epsilon (x)`$ is a normalized, Lorentzian-smoothed $`\delta `$-function of width $`\epsilon `$. The right-hand side of (2) thus corresponds to a sum over eigenstates for which $`E_n`$ lies within a range of size of the order of $`\epsilon `$ centred on $`E`$. Semiclassically, it is approximately the average of $`c_n(𝐋)`$ multiplied by $`\overline{d}(E)`$, the mean level density, which, for systems with two-degrees-of-freedom, is given asymptotically by

$$\overline{d}(E)\frac{V(E)}{(2\pi \mathrm{})^2}$$
(3)
as $`\mathrm{}0`$, where

$$V(E)=\delta (EH(𝐩,𝐪))d^2𝐪d^2𝐩$$
(4)
and $`H(𝐩,𝐪)`$ is the classical Hamiltonian.

Our approach follows the one developed by Keating & Prado (2001) for the special case $`𝐋=\mathrm{𝟎}`$. It is based on the identity

$$G(𝐪+𝐋,𝐪;E)=\underset{n}{}\frac{\psi _n(𝐪+𝐋)\psi _n^{}(𝐪)}{EE_n},$$
(5)
where $`G`$ is the Green function. This implies that

$$\begin{array}{c}\underset{n}{}\psi _n(𝐪+𝐋)\psi _n^{}(𝐪)\delta _\epsilon (EE_n)=\frac{1}{2\pi i}\{G(𝐪+𝐋,𝐪;E+i\epsilon )\hfill \\ \hfill \left(G(𝐪,𝐪+𝐋;E+i\epsilon )\right)^{}\}.\end{array}$$
(6)
For simplicity, we shall concentrate here on systems with time-reversal symmetry. In this case (6) reduces to

$$\underset{n}{}\psi _n(𝐪+𝐋)\psi _n^{}(𝐪)\delta _\epsilon \left(EE_n\right)=\frac{1}{\pi }\mathrm{Im}G(𝐪+𝐋,𝐪;E+i\epsilon ).$$
(7)
(As will become apparent, our analysis generalizes straightforwardly to non-time-reversal-symmetric systems.) The connection with classical mechanics is achieved using the semiclassical approximation to the Green function. For systems with two-degrees-of-freedom, this is

$$G(𝐪+𝐋,𝐪;E)\frac{1}{i\mathrm{}\sqrt{2\pi i\mathrm{}}}\underset{\gamma }{}\sqrt{\left|D_\gamma \right|}\mathrm{exp}\left\{\frac{i}{\mathrm{}}S_\gamma (𝐪+𝐋,𝐪;E)\frac{i\pi }{2}\mu _\gamma \right\},$$
(8)
where the sum includes all classical trajectories from $`𝐪`$ to $`𝐪+𝐋`$ at energy $`E`$, $`S_\gamma `$ is the action along the trajectory labelled $`\gamma `$,

$$D_\gamma =det\left(\begin{array}{cc}\frac{^2S_\gamma }{𝐪^{}𝐪}\hfill & \frac{^2S_\gamma }{𝐪^{}E}\hfill \\ \frac{^2S_\gamma }{E𝐪}\hfill & \frac{^2S_\gamma }{E^2}\hfill \end{array}\right)$$
(9)
with $`𝐪^{}=𝐪+𝐋`$, and $`\mu _\gamma `$ is the Maslov index (Gutzwiller 1990).

The autocorrelation formula we seek follows from averaging (7) with respect to $`𝐪`$. On the left-hand side this gives $`C(𝐋,E,\epsilon )`$. On the right-hand side, the $`𝐪`$-average (via the stationary phase condition) semiclassically selects from the classical orbits with initial and final positions differing by $`𝐋`$ \- henceforth we call such trajectories $`𝐋`$-trajectories - those for which the initial and final momenta are equal. Our main concern will be with cases when $`|𝐋|`$ vanishes in the semiclassical limit. Then the classical orbits contributing to the $`𝐪`$-average are close to periodic orbits. To first approximation one can thus find them by linearizing about the periodic orbits. Essentially, when a periodic orbit is isolated, this corresponds to expanding the action up to terms that are quadratic in the distance from it. Specifically, if $`𝐌_p`$ is the monodromy matrix associated with a periodic orbit $`p`$, then in a section transverse to the orbit, in which the coordinate is $`y`$ and the conjugate momentum is $`p_y`$ (we assume here two degrees of freedom, for simplicity),

$$\left(\begin{array}{c}y+L_y\\ p_y\end{array}\right)=𝐌_p\left(\begin{array}{c}y\\ p_y\end{array}\right),$$
(10)
where $`L_y`$ is the component of $`𝐋`$ in the direction of the $`y`$-axis. For a given choice of $`L_y`$, this allows one to determine the coordinates of the associated orbit in the vicinity of $`p`$ which contributes to $`C(𝐋,E,\epsilon )`$. Thus one may view the autocorrelation function semiclassically as a sum of contributions from the periodic orbits. The contribution from each periodic orbit is, within the linear approximation, stationary in $`L_y`$ when $`p_y=0`$ (i.e. on the periodic orbit itself), that is when $`L_y=0`$. This may be seen explicitly by expanding the action in (7) up to quadratic terms around the orbit $`p`$. Let $`z`$ be the coordinate along $`p`$. Here and in the subsequent analysis we assume that periodic orbit contributions to the autocorrelation function are $`z`$-independent (this is correct in determining the $`\mathrm{}`$ dependence of these contributions when $`|𝐋|0`$ as $`\mathrm{}0`$, because the $`z`$-dependence comes in only via the classical dynamics, except at self-focal points, which do not affect the leading order semiclassical asymptotics of the $`𝐪`$-average). The $`𝐪`$-average may then be computed (it is a Gaussian integral), giving the contribution from $`p`$ to be proportional to $`\mathrm{exp}(i\eta _pL_y^2/\mathrm{})`$, where $`\eta _p`$ is determined by $`𝐌_p`$. It is clear that this is stationary (i.e. centred) at $`L_y=0`$. More importantly for us, it also shows that the contribution from an isolated periodic orbit has a length-scale that is of the order of $`\mathrm{}^{\frac{1}{2}}`$ (in contrast to the random wave model component, for which the length-scale is of the order of $`\mathrm{}`$). This picture coincides precisely with Berry’s (Berry 1989) for Wigner functions. The question we now address is what is the length-scale of the contribution to the autocorrelation function from bifurcating (i.e. non-isolated) periodic orbits.

One might imagine that the autocorrelation formulae for bifurcating orbits could be obtained easily by expanding the action in (8) to higher order than quadratic. This is true for some of the simpler bifurcations (e.g. the codimension-one bifurcations of orbits with repetition numbers $`r=1`$ and $`r=2`$ \- see the examples in Section 3), but for more complicated bifurcations it is incorrect, because for these the linearized map $`𝐌_p^r`$ is equal to the identity, which cannot be generated by the action $`S(𝐪^{},𝐪;E)`$. Thus it is difficult to build into the semiclassical expression for the $`\mathrm{𝐪𝐪}^{}`$-representation of the Green function a well-behaved description of the nonlinear dynamics which the linearized map approximates. The solution to this problem, found by Ozorio de Almeida & Hannay (1987), is to transform the Green function to a mixed position-momentum representation, and this is the approach we now take.

The method just outlined involves, first, Fourier transforming $`G(𝐪^{},𝐪;E)`$ with respect to $`𝐪^{}`$. This gives the Green function in the $`𝐪`$ $`𝐩^{}`$-representation, $`\stackrel{~}{G}(𝐩^{},𝐪;E)`$ ($`𝐩^{}`$ is the momentum conjugate to $`𝐪^{}`$). The semiclassical approximation to $`\stackrel{~}{G}`$ takes the same form as (8), except that $`S(𝐪^{},𝐪;E)`$ is replaced by the $`𝐪`$ $`𝐩^{}`$-generating function $`\stackrel{~}{S}(𝐩^{},𝐪;E)`$. $`G(𝐪^{},𝐪;E)`$ may then be rewritten, semiclassically, as the Fourier transform of this expression with respect to $`𝐩^{}`$. The result, for the semiclassical contribution to $`G(𝐪+𝐋,𝐪;E)`$ from classical orbits in the neighbourhood of a bifurcating periodic orbit, takes the following form.

Consider the general case of a codimension-$`K`$ bifurcation of a periodic orbit with repetition number $`r`$. As before, let $`z`$ be a coordinate along the orbit at bifurcation, let $`y`$ be a coordinate transverse to it, and let $`p_y`$ be the momentum conjugate to $`y`$, so that $`y`$ and $`p_y`$ are local surface of section coordinates. We again assume the semiclassical $`\mathrm{}`$-scaling of the length-scale of the contribution from the bifurcating orbit to the autocorrelation function to be $`z`$-independent. Let $`\mathrm{\Phi }_{r,K}(y,p_y,𝐱)`$ be the normal form which corresponds to the local (reduced) generating function in the neighbourhood of the bifurcation (Arnold 1978; Ozorio de Almeida 1988), where $`𝐱=(x_1,x_2,\mathrm{},x_K)`$ are parameters controlling the unfolding of the bifurcation. Then, up to irrelevant factors, the contribution to $`G(𝐪+𝐋,𝐪;E)`$ is

$$G_{r,K}(y,L_y;𝐱)\frac{1}{\mathrm{}^2}\mathrm{exp}\left[\frac{i}{\mathrm{}}\mathrm{\Phi }_{r,K}(y,p_y,𝐱)+\frac{i}{\mathrm{}}p_yL_y\right]𝑑p_y$$
(11)
(as already stated, we are here interested in determining the $`\mathrm{}`$-dependence of the length-scale of contributions to the autocorrelation function, and so have neglected terms in (11), such as an $`\mathrm{}`$-independent factor in the integrand, which do not influence this). Note that it follows from the definition of the generating function that the stationary phase condition with respect to $`p_y`$ selects classical orbits for which the initial and final $`y`$-coordinates differ by $`L_y`$; that then the stationary phase condition with respect to $`y`$ (which gives the semiclassical contribution to the $`y`$-average needed to compute the autocorrelation formula) further selects from these orbits those whose initial and final momenta, $`p_y`$, are identical; and that the orbits contributing at the point where $`G_{r,K}`$ is stationary with respect to $`L_y`$ are those for which $`p_y=0`$.

We now follow the steps taken in Keating & Prado (2001) for the case when $`L_y=0`$: first, rescale $`y`$ and $`p_y`$ to remove the $`1/\mathrm{}`$ factor from the dominant term (germ) of $`\mathrm{\Phi }_{r,K}`$ in the exponent, and then apply a compensating rescaling of the parameters $`x_1,x_2,\mathrm{},x_K`$ to remove the $`\mathrm{}`$-dependence from the other terms which do not vanish as $`\mathrm{}0`$. Finally, we then rescale $`L_y`$ to remove the $`\mathrm{}`$-dependence of the second term in the exponent in (11). Crucially, this term only involves $`p_y`$, and so if in the initial rescalings

$$p_y\mathrm{}^{\alpha _{r,K}}p_y$$
(12)
then the length scale of the corresponding contribution to the autocorrelation function is of the order of $`\mathrm{}^{1\alpha _{r,K}}`$. This is our main general result. It is important to note that $`\alpha _{r,K}`$ is precisely the exponent found for the scarring amplitude in Keating & Prado (2001). This then implies a sum rule for the scaling exponents of scar amplitudes and the length-scale exponents of wave function correlations reminiscent of those found in the theory of critical phenomena. It is worth emphasizing that the correlation length-scale is related to the amplitude exponent, rather than the width exponent (denoted $`\omega _{r,K}`$ in Keating & Prado 2001).

Values of $`\alpha _{r,K}`$ were computed by Keating & Prado (2001) for the generic bifurcations found in two-degree-of-freedom systems using the appropriate normal forms. For convenience, we list here the results. First, for the $`r=1`$ (i.e. saddle-node) bifurcations which correspond to cuspoid (i.e. corank 1) catastrophes, $`\alpha _{1,K}=\frac{1}{2}`$. Values of $`\alpha `$ for the generic bifurcations with $`r>1`$ and $`K=1`$ are listed in Table 1, and those for $`r>1`$ and $`K=2`$ in Table 2. Finally, $`\alpha _{r,K}=\frac{1}{2(K+1)}`$ in the general case of bifurcations of orbits for which $`r2K+2`$.

Note that in all cases listed above $`0<\alpha 1/2`$, and that in most cases $`\alpha <1/2`$. Thus the length scale of the fringes $`\mathrm{}^{1\alpha _{r,K}}`$ is always semiclassically equal to or smaller than for non-bifurcating orbits, and in most cases is smaller. In all cases it is semiclassically larger than the de Broglie wavelength.

The rescaling of $`y`$ in (11) was defined in Keating & Prado (2001) by

$$y\mathrm{}^{\omega _{r,K}}y,$$
(13)
where $`\omega _{r,K}`$ is the scar width exponent. In the case of the autocorrelation function, one must average over $`y`$. This leads to the amplitude of $`G_{r,K}`$ being proportional to $`\mathrm{}^{\alpha _{r,K}+\omega _{r,K}2}`$. We remark that the exponent of $`\mathrm{}`$ in the amplitude is related to amplitude exponent $`\beta `$ (Berry et al. 2000) of the fluctuations in the spectral counting function associated with the bifurcation in question, because

$$\alpha _{r,K}+\omega _{r,K}2=(\beta _{r,K}+1)$$
(14)
(see equation (2.25) in Keating & Prado 2001).

The rescaling of $`x_n`$ in (11) was defined in Keating & Prado (2001) by

$$x_n\mathrm{}^{\sigma _{n,r,K}}x_n.$$
(15)
Here, the exponents $`\sigma `$ describe the range of influence of the bifurcation in the different unfolding directions $`x_n`$, and their sum

$$\gamma _{r,K}=\underset{n=1}{\overset{K}{}}\sigma _{n,r,K}$$
(16)
describes the $`\mathrm{}`$-scaling of the $`K`$-dimensional $`𝐱`$-space hypervolume affected by the bifurcation.

## 3 Perturbed cat maps

We now illustrate some of the ideas described in the previous section by focusing on a particular example: a family of perturbed cat maps. The influence of bifurcations on the scarring of eigenfunctions was explored numerically for these systems in Keating & Prado (2001). (The influence of tangent bifurcations on localization in wavefunctions was also investigated numerically for a different family of maps by Varga et al. 1999). We focus here on eigenfunction autocorrelations.

The maps we consider are of the form

$$\left(\begin{array}{c}q_{n+1}\hfill \\ p_{n+1}\hfill \end{array}\right)=\left(\begin{array}{cc}2\hfill & 1\hfill \\ 3\hfill & 2\hfill \end{array}\right)\left(\begin{array}{c}q_n\hfill \\ p_n\hfill \end{array}\right)+\frac{\kappa S_p^{^{}}(q_n)}{4\pi ^2}\left(\begin{array}{c}1\hfill \\ 2\hfill \end{array}\right)\text{mod 1},$$
(17)
where $`S_p^{^{}}(q)`$ is the first derivative of a periodic function with period $`1`$, $`\kappa `$ is a parameter determining the size of the perturbation, and $`q`$ and $`p`$ are coordinates on the unit two-torus which are taken to be a position and its conjugate momentum. We will consider two cases: $`S_p(q)=\mathrm{sin}(2\pi q)`$ and $`S_p(q)=\mathrm{cos}(2\pi q)`$.

These maps are Anosov systems. In both cases, for $`\kappa \kappa _{\mathrm{max}}=(\sqrt{3}1)/\sqrt{5}0.33`$ they are completely hyperbolic and their orbits are conjugate to those of the map with $`\kappa =0`$ (i.e. there are no bifurcations). Outside this range, bifurcations occur, stable islands are created, and the dynamics becomes mixed (Berry et al. 1998, Keating & Prado 2001).

The quantization of maps like (17) was developed by Hannay & Berry (1980), when $`\kappa =0`$, and Basílio de Matos & Ozorio de Almeida (1995) for non-zero $`\kappa `$. The quantum kinematics associated with a phase space that has the topology of a two-torus restricts Planck’s constant to taking inverse integer values. The integer in question, $`N`$, is the dimension of the Hilbert space of admissible wavefunctions. With doubly periodic boundary conditions (see, for example, Keating et al. 1999), these wavefunctions in their position representation have support at points $`q=Q/N`$, where $`Q`$ takes integer values between 1 and $`N`$. They may thus be represented by $`N`$-vectors with complex components. The quantum dynamics is then generated by an $`N\times N`$ unitary matrix $`𝐔`$ whose action on the wavefunctions reduces to (17) in the classical limit; for example

$$U_{Q_2,Q_1}=\frac{1}{\sqrt{iN}}\mathrm{exp}\left[\frac{2\pi i}{N}(Q_1^2Q_1Q_2+Q_2^2)+\frac{iN}{2\pi }\kappa \mathrm{sin}(2\pi Q_1/N+\pi \nu /2)\right].$$
(18)
where $`\nu =0`$ gives the sine perturbation and $`\nu =1`$ the cosine perturbation. This matrix plays the role of the Green function of the time-dependent Schrödinger equation for flows.

Denoting the eigenvalues of $`𝐔`$ by $`\mathrm{e}^{i\theta _n}`$ and the corresponding eigenfunctions by $`\mathrm{\Psi }_n(Q)`$, we have, using the fact that the map (17) is invariant under time-reversal and hence that $`𝐔`$ is symmetric, that

$$\begin{array}{c}C(L,\theta ,ϵ)=\underset{n=1}{\overset{N}{}}c_n(L)\delta _\epsilon (\theta \theta _n)=\delta _{L,0}+\frac{1}{N}\mathrm{Re}\underset{Q=1}{\overset{N}{}}\underset{k=1}{\overset{\mathrm{}}{}}U_{Q+L,Q}^k\mathrm{exp}(i\theta k\epsilon k)\hfill \end{array}$$
(19)
where

$$c_n(L)=\frac{1}{N}\underset{Q=1}{\overset{N}{}}\mathrm{\Psi }_n(Q+L)\mathrm{\Psi }_n^{}(Q)$$
(20)
is the eigenvector autocorrelation function and

$$\delta _\epsilon (x)=\frac{1e^\epsilon \mathrm{cos}x}{1+e^{2\epsilon }2e^\epsilon \mathrm{cos}x}$$
(21)
is a periodized, Lorentzian-smoothed $`\delta `$-function of width $`\epsilon `$ (Keating 1991). Equation (19) is the analogue for quantum maps of (7). $`C(L,\theta ,ϵ)`$ corresponds, approximately, to $`N`$ times the local $`n`$-average (over a range of size of the order of $`\epsilon `$) of $`c_n(L)`$. For $`\epsilon `$ large enough, the dominant contribution to (19) comes from the $`k=1`$ term in the sum on the right. We thus define

$$c(L,\theta ,ϵ)=\frac{1}{N}\underset{Q=1}{\overset{N}{}}U_{Q+L,Q}\mathrm{exp}(i\theta \epsilon ).$$
(22)
and, for computational purposes,

$$c_r(L,\theta ,ϵ)=\frac{1}{N}\mathrm{Re}\underset{Q=1}{\overset{N}{}}U_{Q+L,Q}\mathrm{exp}(i\theta \epsilon ).$$
(23)
We may substitute (18) directly into (22). In the semiclassical limit, as $`N\mathrm{}`$, the $`Q`$-average selects regions close to stationary points of the phase of (18). When $`L=0`$ these stationary points coincide with the positions of the fixed points of the classical map (17); otherwise the stationary solutions correspond to $`l`$-trajectories - in this case classical trajectories with initial and final coordinates separated by $`l=L/N`$ and with initial and final momenta equal to each other. These satisfy

$$q_j(l)=\frac{1}{2}\left(jl\frac{\kappa }{2\pi }\mathrm{cos}(2\pi q_j+\nu \frac{\pi }{2})\right)$$
(24)
for integers $`j`$ such that $`0q_j(l)<1`$ (see, for example, Boasman & Keating 1995). When $`l0`$ (which is the most important regime, given that we wish to focus on correlation length-scales that vanish semiclassically) the $`l`$-trajectories obviously lie close to fixed points of the map.

Expanding the phase of (18) around $`q_j`$ up to cubic terms gives

$$\begin{array}{c}U_{Q+L,Q}\frac{1}{\sqrt{iN}}\mathrm{exp}[2\pi iNS_j+\pi iN(2\kappa \mathrm{sin}(2\pi q_j+\nu \frac{\pi }{2}))y^2\hfill \\ \hfill \frac{2\pi ^2iN}{3}\kappa \mathrm{cos}(2\pi q_j+\nu \frac{\pi }{2})y^3],\end{array}$$
(25)
where

$$y=qq_j(l)$$
(26)
and $`2\pi NS_j`$ denotes the phase evaluated at $`q_j(l)`$. Provided that $`2\kappa \mathrm{sin}(2\pi q_j+\nu \pi /2)0`$, this approximation is dominated by the quadratic term in the exponent when $`y`$ is small. It thus describes complex-Gaussian fringes around the classical $`l`$-trajectories with a length-scale (in terms of $`y`$) of the order of $`N^{1/2}`$.

As discussed in Section 2, the centres of the peaks of the autocorrelation function are given by contributing $`l`$-trajectories that start at $`q`$, end at $`q+l`$, and satisfy $`p^{}=p=0`$ in the first iterate of the map (17). From (17), the peak centres are thus solutions of

$$l=q+\frac{\kappa }{2\pi }\mathrm{cos}(2\pi q+\nu \pi /2)n$$
(27)
and

$$3q+\frac{\kappa }{\pi }\mathrm{cos}(2\pi q+\nu \pi /2)m=0$$
(28)
where $`m`$ and $`n`$ are integers. For the linear mapping, when $`\kappa =0`$, the peaks are centred at $`l=q=m/3`$ with $`m=0,1,2`$. These peaks may be seen in Figure 1, where the dashed line represents an evaluation of (23) and the bold line a convolution of the data with a normalized Gaussian of width 0.007.

For $`\kappa >0`$ the peak location equations, (27) and (28), can be solved numerically. As $`\kappa `$ increases bifurcations of periodic orbits take place. It is important at this stage to distinguish between the birth of new periodic orbits as the perturbation parameter $`\kappa `$ is varied (the usual bifurcation) and the birth of new $`l`$-trajectories as $`l`$ varies (a process which we call an $`l`$-caustic). As $`l0`$, the $`l`$-trajectories are close to periodic orbits. $`l`$-caustics coincide with periodic orbit bifurcations in the limit.

Our main purpose now is to illustrate quantitatively the influence of bifurcating periodic orbits on the autocorrelation function. We shall study in detail two bifurcations: one a first order (tangent) bifurcation, and the other a second order bifurcation. In order to do so, we will treat separately the sine and cosine perturbations.

### 3.1 Sine perturbation $`\nu =0`$

For $`\kappa >0`$, the positions of the peaks in the autocorrelation function deviate from those ($`0`$, $`1/3`$ and $`2/3`$) of the unperturbed map ($`\kappa =0`$). The real part of the autocorrelation function for $`\kappa =1`$ is shown in Figure 2 (a). The peak close to $`l=0`$ comes from the contribution of the period-1 fixed point $`j=0`$ in (24). Using a stationary phase approximation with respect to $`L`$ in (25) gives the position of the peak to be

$$l\frac{\kappa }{\pi }\frac{\mathrm{cos}(2\pi q_0(0))}{3\kappa \mathrm{sin}(2\pi q_0(0))}\text{ mod 1}$$
(29)
in the limit as $`l0`$. For $`\kappa \mathrm{sin}(2\pi q_j))<2`$, the second order term in the Taylor expansion (25) dominates and the semiclassical formula for the autocorrelation function is simply

$$c_r(L,\theta ,ϵ)\underset{j=0}{\overset{1}{}}\mathrm{exp}(\epsilon )\frac{\mathrm{cos}(2\pi NS_j\theta +(\zeta 1)\pi /4)}{\sqrt{|2\kappa \mathrm{sin}(2\pi q_j(l))|}}$$
(30)
where $`\zeta =\text{sign}(2\kappa \mathrm{sin}(2\pi q_j(l)))`$ and

$$S_j=\frac{3l^2}{4}+j\frac{l}{2}j\frac{1}{4}+\frac{\kappa }{4\pi ^2}\mathrm{sin}(2\pi q_j(l))+\frac{\kappa ^2}{16\pi ^2}\mathrm{cos}^2(2\pi q_j(l)).$$
(31)
In figure 2 (b) we plot the difference between the exact expression (23) and the semiclassical approximation (30).

For larger $`\kappa `$, keeping terms only up to second order in the approximation (25) fails in two ways. First, $`2\kappa \mathrm{sin}(2\pi q_j(l))`$ vanishes for some values of $`l`$ and higher terms of the expansion of the action need to be taken into account. It is clear in Figure 3, where the predicted peak-centres (obtained by solving (27) and (28)) are plotted, that for $`\kappa `$ around $`3.0`$ new $`l`$-trajectories contribute to the autocorrelation with a peak around $`l=0.5`$. Although these $`l`$-caustics are an additional complication, they are in fact unimportant if the object is to compute the contribution of bifurcating periodic orbits which takes place for $`l`$ around zero (recall that this is the region that is semiclassically interesting for correlation length-scales that vanish with $`\mathrm{}`$).

Second, when $`l0`$, $`2\kappa \mathrm{sin}(2\pi q_j(0))=0`$ at periodic orbit bifurcations. Then again the quadratic term in (25) vanishes, and the fringe structure comes from the cubic term. It thus has a $`y`$-length-scale of the order of $`N^{1/3}`$. In the language of Section 2, this corresponds to a codimension-one bifurcation of a periodic orbit with $`r=1`$ (a tangent bifurcation). The $`j=0`$ and $`j=1`$ terms each correspond to a single unstable fixed point if the condition $`2\kappa \mathrm{sin}(2\pi q_j(l))0`$ is satisfied. Note that at both $`l`$-caustics and tangent bifucations, the semiclassical asymptotics is determined by the cubic term in the expansion of the action; hence both contributions are of the same order in $`\mathrm{}`$.

Figure 4 provides a panorama of the situation. Here, the absolute value of the Gaussian smoothed autocorrelation function $`c`$ is plotted. The result should be compared with the predicted peak-centres shown in Figure 3. The peak coming from the bifurcating periodic orbits is around $`\kappa =6`$ and $`l=0`$. The other similar peaks are all associated with the $`l`$-caustics. It is clear that their amplitudes are much larger than those of the peaks associated with isolated $`l`$-trajectories.

For the sine perturbation, the first bifurcation occurs when $`\kappa =\kappa ^{}=5.943388`$. At this parameter value, two new degenerate solutions of (24) appear for both $`j=0`$ and $`j=1`$, corresponding in each case to the birth of a pair of fixed points, one stable and the other unstable. For $`l=0`$ the contributions of these bifurcating orbits when $`\kappa =\kappa ^{}`$ are:

$$c_r(0,\theta ,ϵ)\underset{j=0}{\overset{1}{}}\frac{N^{1/6}}{\mathrm{\Gamma }(2/3)}\left(\frac{4\pi }{9\kappa ^{}}\right)^{1/3}\mathrm{exp}(\epsilon )\frac{\mathrm{cos}(2\pi iN\overline{S}_j\theta \pi /4)}{|\mathrm{cos}(2\pi \overline{q}_j)|^{1/3}}$$
(32)
where $`\overline{S}_j=(S(q_{+}^{}{}_{j}{}^{})+S(q_{}^{}{}_{j}{}^{}))/2`$ is the average of the actions of the two bifurcating fixed points at $`q_{+}^{}{}_{j}{}^{}`$ and $`q_{}^{}{}_{j}{}^{}`$. More generally, the semiclassical contributions of the real orbits plus the complex orbits to the autocorrelation function give:

$$\begin{array}{c}c_r(L,\theta ,ϵ)\underset{j=0}{\overset{1}{}}\mathrm{exp}(\epsilon )\frac{\mathrm{cos}(2\pi NS_j\theta +(\sigma 1)\pi /4)}{\sqrt{|2\kappa \mathrm{sin}(2\pi q_j(l))|}}+\hfill \\ \hfill \mathrm{Re}\underset{j=0}{\overset{1}{}}\mathrm{exp}(\epsilon )N^{1/6}\frac{\text{Ai}\left[(2\pi N|\alpha |)^{2/3}\right]}{2\pi |\alpha |^{1/3}}\sigma ^2\mathrm{exp}(2\pi iN\overline{S}_ji\theta i\pi /4).\end{array}$$
(33)
Here $`\sigma =(q_+q_{})/2`$ and $`|\alpha |=3|\mathrm{\Delta }S|/(2\sigma ^3)`$, where $`\mathrm{\Delta }S=(S(q_+)S(q_{}))/2`$. Figure 5 (a) shows the autocorrelation function for $`\kappa =\kappa ^{}`$. In Figure 5 (b) the deviation of the semiclassical approximation (33) from the exact expression (23) is plotted.

### 3.2 Cosine perturbation $`\nu =1`$

Setting $`\nu =1`$ in (18) and the following equations allows us now to discuss a second order bifurcation. (Such bifurcations occur for the second iterate of the map with the sine perturbation, and so quantum mechanically influence $`U^2(Q,Q)`$. In the case of the cosine perturbation, they occur for fixed points of the map, and so influence $`U(Q,Q)`$. This makes them easier to study numerically (Mende (1999).)

As the order of the first non-vanishing derivative of the action $`S_p(q)`$ determines the order of the bifurcation, one can see that there are no higher bifurcations than those of second order for this perturbation. For $`j=0`$, $`q=0`$ and $`\kappa =2`$, a bifurcation of order two occurs. For $`j=1`$ there is no such bifurcation: in this case only tangent bifurcations occur, the first one at $`\kappa 9.208`$ (Mende 1999). As the tangent bifurcation was discussed previously, we will concentrate on the bifurcation of order two. For this, when $`\kappa <2`$ there is one real unstable fixed point (called the central fixed point), while for $`\kappa >2`$ there are two real fixed points (satellites) in addition to the central fixed point, which is now stable. The scarring of the wavefunctions caused by this bifurcation is shown as a dashed line is Figure 6, where $`_n|\mathrm{\Psi }_n(Q)|^2\delta _\epsilon (\theta \theta _n)1`$ is plotted for (a) $`\kappa =1`$ (one real unstable fixed point ($`j=0`$) at $`q=0`$ and another real unstable fixed point ($`j=1`$) at $`q=0.5`$), (b) $`\kappa =\kappa ^{}=2`$ (where the $`j=0`$ fixed point bifurcates), and (c) $`\kappa =3`$ (three real $`j=0`$ fixed points - unstable satellites at $`q=0.23`$ and $`q=0.23`$ and the central fixed point at $`q=0`$, which is now stable - and the unstable $`j=1`$ fixed point at $`q=0.5`$). A convolution of the data with a normalized Gaussian of width 0.02 is also shown (bold line).

In Figure 7, $`|_{Q=1}^NU(Q,Q)|^2=|\mathrm{Tr}U|^2`$ is plotted as a function of $`\kappa `$ and $`N`$. One sees the contribution of the second order bifurcation of the $`j=0`$ fixed point around $`\kappa =2`$ and that of the tangent bifurcation of the $`j=1`$ fixed point around $`\kappa =9`$. The amplitude associated with the second order bifurcation is clearly semiclassically larger. Similarly, the amplitude of the autocorrelation function is larger at these points (compare (34) below with (32)). The predicted positions of the peaks in the autocorrelation function (obtained by solving (27) and (28), when $`\nu =1`$) are shown in Figure 8. Here again, for $`\kappa =0`$ the peaks are centred at $`l=01`$, $`l=1/3`$ and $`l=2/3`$, as previously described. As $`\kappa `$ increases, these peaks deviate from their original positions, but the autocorrelation result for isolated points, the analogue of (30), is still valid. For $`\kappa `$ around 2, where the second order bifurcation occurs, the approximation (25) breaks down and an expansion of the phase of (18) around $`q_j`$ up to fourth order terms is necessary. The contribution of the bifurcating fixed point ($`j=0`$) when $`\kappa =\kappa ^{}=2`$ and $`L=0`$ is

$$c_r(0,\theta ,ϵ)\frac{N^{1/4}}{2\sqrt{2}\mathrm{\Gamma }(3/4)}\left(\frac{3\pi }{\kappa ^{}}\right)^{1/4}\mathrm{exp}(\epsilon )\mathrm{cos}(2\pi iNS_0\theta \pi /4)$$
(34)
where $`S_0`$ is the action of central orbit at $`q=0`$. More generally, for $`\kappa 9`$ the semiclassical approximation to the autocorrelation function is

$$\begin{array}{c}c_r(L,\theta ,ϵ)\{\frac{\mathrm{cos}(2\pi NS_1\theta +(\zeta 1)\pi /4)}{\sqrt{|2\kappa \mathrm{cos}(2\pi q_1(l))|}}\hfill \\ \hfill +\mathrm{Re}\frac{\pi }{2}N^{1/2}q_s(l)[\text{J}_{1/4}(z)\mathrm{exp}(i\sigma _1)+\text{sign}(2\kappa )\text{J}_{1/4}(z)\mathrm{exp}(i\sigma _2)]\}\mathrm{exp}(\epsilon ).\end{array}$$
(35)
Here $`q_s(l)`$ is the distance of the satellite points from the origin, $`z=\pi N\alpha q_s^4/4`$, where $`\alpha =4(q_s^2+\kappa (1\mathrm{cos}(2\pi q_s)))/q_s^4`$, $`\sigma _1=\theta \pi /8+2\pi NS_0\pi N\alpha q_s^4/4`$, and $`\sigma _2=\theta +5\pi /8+2\pi NS_0\pi N\alpha q_s^4/4`$.

In Figure 9 (a) the real part of the autocorrelation function for $`\kappa =\kappa ^{}=2`$ is plotted and in (b) the deviation of the semiclassical approximation (35) from the exact expression (23) is shown. In Figure 10 the absolute value of the convolution of the autocorrelation function (22) with a normalized Gaussian is plotted as a function of $`\kappa `$ and $`l`$, showing the consistency of the theory: higher amplitudes near $`l`$-caustics and closed agreement with the predicted positions of the peak-centres (Figure 8). Finally, in Figure 11 the absolute value of the real part of the autocorrelation function (23) at $`L=0`$ (convolved with a Gaussian in $`\kappa `$ of width 0.1) is plotted. From this one sees clearly the influence of the two bifurcations: the tangent bifurcation at $`\kappa 9.208`$ and the second order bifurcation at $`\kappa =2`$.

## 4 Acknowledgements

AB thanks the School of Mathematics, University of Bristol, and the Basic Research Institute in the Mathematical Sciences (BRIMS), Hewlett-Packard Laboratories, Bristol, for hospitality during the period when this work was started. SDP is grateful to Martin Sczyrba for helpful discussions, and to the School of Mathematics at the University of Bristol for hospitality and financial support. Support is also acknowledged from the Deutsche Forschungsgemeinschaft (for AB) and the Fundação de Amparo à Pesquisa do Estado do Rio Grande do Sul - FAPERGS (for SDP).

References

Arnold, V.I. 1978 Mathematical Methods in Classical Mechanics. Springer.

Bäcker, A. & Schubert, R. 2002 Amplitude distribution of eigenfunctions in mixed systems. J. Phys. A 35, 527-538 (BS1).

Bäcker, A. & Schubert, R. 2002 Autocorrelation function of eigenstates in chaotic and mixed systems. J. Phys. A 35, 539-564 (BS2).

Basílio de Matos, M. & Ozorio de Almeida, A.M. 1995 Quantization of Anosov maps. Ann. Phys. 237, 46-65.

Berry, M.V. 1977 Regular and irregular semiclassical wavefunctions J. Phys. A 10, 2083-2091.

Berry, M.V. 1989, Quantum scars of classical closed orbits in phase space. Proc. R. Soc. Lond. A 243, 219-231.

Berry, M.V, Keating, J.P & Prado, S.D. 1998 Orbit bifurcations and spectral statistics. J. Phys. A 31, L245-254.

Berry, M.V., Keating, J.P. & Schomerus, H. 2000 Universal twinkling exponents for spectral fluctuations associated with mixed chaology. Proc. R. Soc. Lond. A 456, 1659-1668.

Boasman, P.A. & Keating, J.P. 1995 Semiclassical asymptotics of perturbed cat maps. Proc. R. Soc. Lond. A 449, 629-653.

Bogomolny, E.B. 1988 Smoothed wavefunctions of chaotic quantum systems. Physica D 31, 169-189.

Gutzwiller, M.C. 1990 Chaos in Classical and Quantum Mechanics (New York: Springer).

Hannay, J.H. & Berry, M.V. 1980 Quantization of linear maps on the torus – Fresnel diffraction by a periodic grating. Physica D 1, 267-290.

Heller, E.J. 1984 Bound state eigenfunctions of classically chaotic Hamiltonian systems - scars of periodic orbits. Phys. Rev. Lett. 53, 1515-1518.

Kaplan, L. 1999 Scars in quantum chaotic wavefunctions. Nonlinearity 12, R1-R40.

Keating, J.P. 1991 The cat maps: quantum mechanics and classical motion. Nonlinearity 4, 309-341.

Keating, J.P., Mezzadri, F. & Robbins, J.M. 1999 Quantum boundary conditions for torus maps. Nonlinearity 12, 579-591.

Keating, J.P. & Prado, S.D. 2001 Orbit bifurcations and the scarring of wave functions. Proc. R. Soc. Lond. A 457, 1855-1872.

Li, B. & Rouben, D.C. 2001 Correlations of chaotic eigenfunctions: a semiclassical analysis. J. Phys. A 34, 7381-7391.

Mende, M. 1999 Periodic orbit bifurcations and spectral statistics. M.Sc. Thesis, University of Bristol.

Ozorio de Almeida, A.M. 1988 Hamiltonian Systems: Chaos and Quantization. Cambridge University Press.

Ozorio de Almeida, A.M. & Hannay, J.H. 1987 Resonant periodic orbits and the semiclassical energy spectrum. J. Phys A 20, 5873-5883.

Percival, I.C. 1973 Regular and irregular spectra. J. Phys. B 6, L229-L232.

Varga, I., Pollner, P. & Eckhardt, B. 1999 Quantum localization near bifurcations in classically chaotic systems. Ann. Phys. (Leipzig) 8, SI265-SI268.


format: text

text:
# Extended Star Formation in Dwarf Spheroidal Galaxies: the Cases of Draco, Sextans, and Ursa Minor

## 1 Introduction

Satellite dwarf spheroidal galaxies of the Milky Way have been considered to be similar to globular clusters (e.g., Hodge 1971), since they are composed of old metal-poor stars and contain little gas. Recent analyses of the colour-magnitude diagrams (CMDs), however, reveal a surprising fact that dwarf spheroidal galaxies (dSphs) are of composite stellar populations. Some of dSphs formed stars continuously and/or intermittently for many Gyr (e.g., the dSph Carina: Smecker-Hane et al. 1996; the dSph Leo I: Gallart et al. 1999). Low-resolution spectra of individual stars along the red giant branch (RGB) have shown that both Galactic and M31 satellite dSphs have internal metallicity dispersions (e.g., Suntzeff et al. 1993; Smecker-Hane & McWilliam 1999; C$`\widehat{\mathrm{o}}`$t$`\stackrel{´}{\mathrm{e}}`$, Oke, & Cohen 1999). Using high dispersion (R=34000) spectroscopy, Shetrone, C$`\widehat{\mathrm{o}}`$t$`\stackrel{´}{\mathrm{e}}`$, & Sargent (2001) have recently presented abundances of various elements in the three dSphs Draco, Sextans, and Ursa Minor, which are all satellite galaxies of the Milky Way. They report that the abundance ratios of $`\alpha `$ elements (O, Si, Mg, etc.) relative to iron (\[$`\alpha `$/Fe\]) are nearly equal to solar at low metallicity (\[Fe/H\] $`<1.4`$). This suggests that the stars were formed in the gas enriched by Type Ia supernovae (SNe Ia) as well as by Type II supernovae (SNe II), that is, star formation in the dSphs continued over the lifetime (typically $`12`$ Gyr) of SNe Ia progenitors. This result allows us to directly study the chemical enrichment histories (CEHs) of the dSphs Draco, Sextans, and Ursa Minor.

One problem in deriving an accurate SFH from stellar colours is the so-called age-metallicity degeneracy. Older stars become redder and the same is true for more metal-rich stars. Hence, young and metal-rich stars show colours similar to old and metal-poor ones. We have empirically solved the age-metallicity degeneracy for the three dSphs cited above to derive their star formation histories (SFHs) to a level of accuracy which is impossible in other galaxies in the Local Group. To reveal SFHs in detail would help to understand how large galaxies and the environments affect the evolution of dwarf galaxies.

In addition to a study of the evolution of dwarf galaxies themselves, one can also learn about the possible merging history of the Milky Way by comparing the abundances of the satellite dwarf galaxies with those of Galactic field halo stars. In a hierarchical galaxy formation scenario, larger galaxies grew at the expense of smaller gaseous fragments (e.g., White & Rees 1978; Kauffmann, White & Guiderdoni 1993; Blumenthal et al. 1994; Cole et al. 1994). From the observational perspective, Searle & Zinn (1978) proposed a model in which the Galactic halo was formed via infall and destruction of proto-Galactic fragments. Increasing empirical and analytical evidence suggests that the Galactic halo was, at least in part, assembled from chemically-distinct, low-mass fragments (e.g., Searle & Zinn 1978; Yanny et al. 2000). Such evidence includes the recent discovery of the tidally disturbed dSph Sagittarius which is falling onto the Galaxy (Ibata, Gilmore & Irwin 1994) as well as the numerous reports of kinematical substructure among halo field stars (e.g., Chen 1998). The scenario can be tested by comparing abundance patterns of field stars in the Galactic halo and those of satellite dwarf galaxies. If today’s dwarf galaxies in the Local Group were the counterparts of proto-Galactic gaseous fragments, the abundance patterns of dwarf galaxies should be similar to those of field stars in the Galactic halo. The observations by Shetrone et al. (2001), however, have shown that the patterns of abundance ratio in the dSphs are different from those of the Galactic halo stars. The dSphs have $`0.4[\mathrm{Mg}/\mathrm{Fe}]0.4`$, while field stars in the Galactic halo have $`0.2[\mathrm{Mg}/\mathrm{Fe}]0.6`$ . The authors concluded that the Galactic halo was unlikely to have formed via the accretion of objects similar to the dSphs Draco, Sextans, and Ursa Minor.

In this paper, we discuss the SFHs and the CEHs of the three dSphs Draco, Sextans, and Ursa Minor. Using the abundance patterns and the CMDs, we show that the low (1$`5`$ % of the solar neighbourhood) star formation rate and relatively long duration ($`3.96.5`$ Gyr) of initial star formation explain both the abundance patterns and morphologies of CMDs such as well-populated red horizontal branch (RHB) and the tight red giant branch (RGB).

The paper is organised as follows. In section 2, chemical evolution model is confronted with the observed stellar abundances to study the CEHs and SFHs. In section 3, by calculating CMDs by using our CMD simulation code (Ikuta 2001), we will show that the SFHs derived in the previous section succeed in reproducing observed CMDs. Discussions and conclusions are presented in sections 4 and 5, respectively.

## 2 Observed Abundances and Chemical Evolution

### 2.1 Chemical evolution model

Generally stellar birth rate is separated into two independent functions. The birth rate of stars with mass between $`m`$ and $`m+dm`$ is described as $`C(t)\varphi (m)dm`$, where $`C(t)`$ and $`\varphi (m)`$ are the SFR and the initial mass function (IMF), respectively (Tinsley 1980). The IMF is assumed to be time invariant with a power-law spectrum. Normalising to unity, we have

$$\varphi (m)=\frac{(x1)m_l^{x1}}{1(m_l/m_u)^{x1}}m^x,$$
(1)
where the lower and the upper mass limits are assumed to be $`m_l=0.1M_{}`$ and $`m_u=60M_{}`$, respectively. The Salpeter IMF has $`x=1.35`$ in this definition. Similarly to the solar neighbourhood disc model of chemical evolution (Arimoto, Yoshii, & Takahara 1992), the Salpeter IMF is assumed in this section. Models with different IMFs will be discussed in section 4. The gas mass $`M_\mathrm{g}(t)`$ changes through star formation:

$$\frac{dM_\mathrm{g}(t)}{dt}=C(t)+E(t),$$
(2)
where $`E(t)`$ is the gas ejection rate from dying stars, assuming a closed-box model. In the case of the Galactic disc, an infall model is generally assumed to solve the so-called G-dwarf problem (e.g., Tinsley 1980; Arimoto et al. 1992). However, there is no observational evidence for the G-dwarf problem in dwarf galaxies, thus we do not consider the infall of gas here. Gas outflow during the star formation is not considered, either. A recent numerical experiment showed that the fraction of mass loss from dwarf galaxies is small (at most $`7`$ percent of the galaxy mass; MacLow & Ferrara 1999). Metal-enhanced wind is often introduced to explain observed properties of dwarf galaxies (e.g., Vader 1986). MacLow & Ferrara (1999) showed that indeed most of the metals are able to leave dwarf galaxies during star burst, although the mass loss rate is small. Virtually no metals are retained in dwarf galaxies if high energy is supplied by supernova explosions (see table 3 in their literature for details). Complete metal loss means no chemical enrichment in dwarf galaxies. This is inconsistent with observations which show clearly that dSphs are chemically polluted. Their results suggest that only star formation of a low rate allows a dwarf galaxy to retain a large fraction of heavy elements produced in massive stars. As will be shown later, very low SFRs are indeed required to explain the observed abundances and CMDs. Thus, a simple model is adopted here and a galaxy is assumed to be a closed system.

The SFR is assumed to be proportional to the gas fraction $`f_\mathrm{g}(t)`$:

$`C(t)=ϵ_{\mathrm{SF}}M_\mathrm{G}f_\mathrm{g}(t)=ϵ_{\mathrm{SF}}M_\mathrm{g}(t)`$ (3)
with

$$f_\mathrm{g}(t)M_\mathrm{g}(t)/M_\mathrm{G},$$
(4)
where $`ϵ_{\mathrm{SF}}^1`$ (Gyr) and $`M_\mathrm{G}`$ are the timescale of star formation and the initial total mass of a galaxy, respectively.

The evolution of the abundance of the $`i`$-th element $`Z_i(t)`$ is given as

$$\frac{d(Z_iM_\mathrm{g})}{dt}=Z_iC(t)+E_{Z_i}(t)+Z_i(0)M_\mathrm{G},$$
(5)
where $`E_{Z_i}(t)`$ is the total ejection rate of processed and unprocessed $`i`$-th element, $`Z_i(0)`$ is the initial abundance of $`i`$-th element in a proto-cloud. Hereafter, $`Z_i(0)=0`$ is always assumed, i.e., a proto-cloud was metal-free. We calculate $`E(t)`$ and $`E_{Z_i}(t)`$ from the following integrals:

$$E(t)=_{m_t}^{m_u}(1w_m)C(tt_m)\varphi (m)𝑑m,$$
(6)
$`E_{Z_i}(t)`$ (7)
$`={\displaystyle _{m_t}^{m_u}}\left[(1w_m)Z_i(tt_m)+p_{im}\right]C(tt_m)\varphi (m)𝑑m,`$
where $`t_m`$ is the lifetime of a star with mass $`m`$ and the lower limit $`m_t`$ is the stellar mass with lifetime $`t_m=t`$. Nucleosynthesis data $`w_m`$ and $`p_{im}`$ are the remnant mass fraction and the total mass fraction of processed and unprocessed $`i`$-th element. It is still debated whether stellar metallicity affects both $`w_m`$ and $`p_{im}`$, but no firm consensus has been reached yet. Thus we neglect these possible effects: instead we adopt the values calculated for the solar abundance. We assume that the lower mass limit $`m_{l,\mathrm{II}}`$ for a progenitor of Type II supernova (SNeII) is $`10M_{}`$ (Tsujimoto et al. 1995). Nucleosynthesis of SNeII and of SNeIa are explicitly included, since the former profoundly produce the $`\alpha `$ elements (e.g., O, Mg, Si, Ca), while the latter the iron-peak elements (e.g., Fe, Ni, Co). The data of nucleosynthesis are taken from Thielemann, Nomoto, & Hashimoto (1996) for SNeII, Nomoto, Thielemann, & Yokoi (1996) for SNeIa, and Renzini & Voli (1981) for low-mass ($`8M_{}`$) stars. Progenitors of SNe Ia are generally considered to have longer lifetimes than those of SNe II, since a majority of them are likely to be Chandrasekhar mass white dwarfs (e.g., Nomoto, Iwamoto, & Kishimoto 1997). Observations for disc (Edvardsson et al. 1993) and halo (e.g., McWilliams et al. 1995; Gratton & Sneden 1987) stars in the Milky Way show that there is a sudden decline of \[$`\alpha `$/Fe\] in the \[$`\alpha `$/Fe\]-\[Fe/H\] diagram. The decline is generally interpreted as an onset of SN Ia explosions. Studying observed \[O/Fe\] break at \[Fe/H\] $`1`$ for disc and halo stars in the Milky Way, Yoshii, Tsujimoto, & Nomoto (1996) suggested that the effective lifetime of SNe Ia is around $`1.5`$ Gyr. A new SNIa model (Hachisu, Kato, & Nomoto 1999) suggested that occurrence of SNeIa depends on the metallicities of progenitors as well as the lifetimes. Based on the new SN Ia model, Kobayashi et al. (1998) predicted that SN Ia events occur only for the progenitors with \[Fe/H\] $`1.1`$. On the observational side, $`\alpha `$ element to iron abundance ratios of the three dSphs, Draco, Sextans, and Ursa Minor (Shetrone et al. 2001) have indicated a significant contribution from SNeIa at a metallicity lower than \[Fe/H\]$`=1.0`$. The observed trends of \[$`\alpha `$/Fe\] are difficult to explain by the new SNIa model. Thus, the metallicity effect on SNeIa is not considered here and the lifetime of the progenitors is assumed to be $`1.5`$ Gyr (Yoshii et al. 1996).

### 2.2 Comparison with the observation

Fig. 1 reproduces the observed stellar abundances of the dSphs Draco, Sextans, and Ursa Minor in the \[Mg/Fe\]$``$\[Fe/H\] diagram (Shetrone et al. 2001). Despite the narrow widths of the RGB (see figure 1 of Shetrone et al. 2001), the iron abundances exhibit a large dispersion ($`3[\mathrm{Fe}/\mathrm{H}]1.4`$), and the abundance ratios of the $`\alpha `$-elements (e.g., Mg) to iron are near or below solar. The abundance ratio \[Mg/Fe\] starts to decrease at \[Fe/H\] $`2`$. On the other hand, \[Mg/Fe\] in the solar neighbourhood decreases with \[Fe/H\] at around \[Fe/H\]$`1`$, which is generally interpreted as the onset of SNeIa explosions (e.g., Matteucci & Greggio 1986). Supposing that the decline of \[Mg/Fe\] in the dSphs is also caused by SNeIa, the star formation should have continued longer than the lifetime (typically $`12`$ Gyr) of progenitors of SNeIa and the SFRs should be much lower than that of the solar neighbourhood. A time-delay model (Matteucci & Brocato 1990) of the SNeIa precisely predicted that the low star formation rate results in lower \[Mg/Fe\] ratios relative to stars in the Galaxy.

In Fig. 1, theoretical abundance patterns of our models A$``$D are superposed. The models assume different $`ϵ_{\mathrm{SF}}`$ of $`210^2,110^2,210^3,110^3`$ (Gyr-1) as indicated in Fig. 1. The model parameters are summarised in table 1. A typical value of $`ϵ_{\mathrm{SF}}`$ for the solar neighbourhood disc is $`ϵ_{\mathrm{SF}}=0.2`$ Gyr-1 (Arimoto et al. 1992).

The models B$``$D shown in Fig. 1 give a good fit to the observed abundance patterns. Since the abundance ratios resulting from the model A are too high to reproduce the observed ones, the model A is rejected. The model D with $`ϵ_{\mathrm{SF}}=110^3`$ Gyr-1 is not appropriate, either. Since chemical evolution is so slow, it takes around 15 Gyr to reach \[Fe/H\]$`=1.4`$, which is the highest abundance observed in these galaxies. If this model is adopted, the star formation must continue till now. Clearly this is inconsistent with the fact that no star formation occurs currently in these galaxies. Thus, only models B and C remain as being acceptable. We note that the observational data for the dSph Sextans have lower signal-to-noise ratios than those for the dSphs Draco and Ursa Minor (see table 2 of Shetrone et al. 2001) and the abundances of the dSph Sextans are less reliable than those of the other two dSphs. Thus, it is not necessary to take the discrepancy between the model predictions and the abundances of the dSph Sextans too seriously.

In short, the abundance patterns of the dSphs Draco, Sextans, and Ursa Minor all suggest that the chemical enrichment occurred with SFRs which are much lower than normal spiral galaxies (only $`15\%`$).

## 3 CMDs

The models B and C fit well to the observed trends in the \[Mg/Fe\]$``$\[Fe/H\] diagram. Next, we discuss if the colour-magnitude diagrams (CMDs) predicted by these models can reproduce the observed ones. The CMDs provide us with the most detailed information to follow SFHs and CEHs back to the oldest stars. The SFH derived from the abundance patterns in the previous section must give a consistent CMD with the observations.

### 3.1 Simulation code of a CMD

Simulating numerical CMDs has become a standard technique to study the SFHs of nearby galaxies through the observed CMDs. Galaxies are composed of complex stellar populations, and the observed CMDs are affected by photometric errors and detection incompleteness. A Monte Carlo simulation allows us to take into account both simultaneously. A composite stellar population is randomly generated according to predictions of stellar evolution theory with the IMF and SFR being assumed, and then uncertainties of the data such as the increasing scatters and rising incompleteness at fainter magnitudes are taken into account.

Because of these advantages, this approach is becoming more and more frequently used to study the stellar populations of nearby galaxies. The approach was first applied by Ferraro et al. (1989) and was fully described by Tosi (1991) and Greggio et al. (1993). A more quantitative method, the R-method, to compare simulated and observed CMDs was presented by Bertelli et al. (1992) and was described by Vallenari et al. (1996) in more detail. The statistical comparison was further developed by Tolstoy & Saha (1996) and Tolstoy (1996) by adopting Bayes’ theorem. Now several groups have constructed CMD simulators and have investigated SFHs of galaxies in the Local Group (e.g., Aparicio et al. 1996; Dolphin 1997; Hernandez, Valls-Gabaud, & Gilmore 1999; Gallart et al. 1999).

We should point out, however, that there were two problems in the previous studies. First, all the approaches so far developed use the so-called optimising method of stellar population synthesis. The best mixture of stellar population is searched iteratively to reproduce the CMDs. Problems with this approach are the unproven uniqueness of the solution (e.g., Greggio et al. 1999) and a lack of evolutionary information.

Second, stellar metallicity was assumed a priori in all the simulation codes presented previously. Some introduced metallicity variation in time (e.g., Gallart et al. 1999), but the variation assumed was independent of the SFH. Since the metallicity affects evolutionary tracks and atmospheres of stars, the colours and luminosities of stars are changed by metallicity. Thus, it is crucial to use a simulator of CMDs which fully takes into account a CEH to interpret the CMDs properly and to derive the SFHs accurately. The simplification and/or neglect of chemical evolution lead to a serious problem called the age-metallicity degeneracy. Since stellar colours become bluer when stars are younger and/or poorer in metallicity, there are at least two interpretations for a given position of a star in the CMD; young and metal-rich or old and metal-poor. Therefore, an assumption of the metallicity (or age) of a certain stellar population may result in wrong estimation of age or metallicity, due to the age-metallicity degeneracy.

Aiming to solve the degeneracy, we have built a numerical simulation code of CMD morphology (Ikuta 2001). In the code, we adopt an evolutionary method of stellar population synthesis (e.g., Arimoto & Yoshii 1986). A galaxy is assumed to have been a proto-galactic gas cloud at the beginning. Stars formed and newly processed elements were ejected from stars at the end of their lives either via stellar wind or via supernova explosions (SNeIa and SNeII). The gas was chemically enriched and the next generation of stars changed their evolution and photometric properties due to a metallicity increase. The code is particularly designed to disentangle stellar age and metallicity both of which heavily affect the morphology of CMDs. To take into account effects of chemical evolution, a fine interpolation in stellar ages and metallicities is adopted in our code. The later stages of stellar evolution, in particular, horizontal branch (HB), asymptotic giant branch (AGB), and post-AGB stars are treated in a sophisticated way, since these bright stars are crucial for this purpose. In this section, the procedure of our Monte Carlo simulation of CMDs and the ingredients are detailed.

#### 3.1.1 Outline of the Simulation Procedure

Fig. 2 shows the flow-chart of the simulation. First, a SFH is assumed. The parameters to describe the SFH are the timescale of star formation $`ϵ_{\mathrm{SF}}^1`$ (Gyr), the duration of star formation $`\mathrm{\Delta }T_{\mathrm{SF}}`$ (Gyr), the slope of the IMF ($`x`$), and the age of a galaxy. Second, chemical evolution is calculated, which gives the total number of stars at each metallicity at each time step. Third, a HR diagram is calculated by using stellar evolutionary tracks, where a Monte Carlo simulation is performed to assign stellar masses. Fourth, the effective temperature $`T_{\mathrm{eff}}`$ and luminosity $`L`$ of a star are transformed to colour and magnitude by using a library of stellar model atmospheres. Fifth, a Monte Carlo simulation is performed again to take into account observational errors and detection completeness, and thus the CMD is obtained. A number of stars in the simulated CMD is adjusted to an observed CMD. Finally, the simulated CMD is confronted with the observed one by visual inspection. If they are inconsistent, the parameters of the SFH are changed, and the simulation is re-started. If the simulated CMD fits the observation, the SFH and CEH are derived. Ingredients of the simulation, such as stellar evolutionary tracks, model atmospheres, and observational conditions, are described in the subsequent section.

#### 3.1.2 Stellar Tracks and Model Stellar Atmosphere

Stellar evolutionary tracks give luminosities, effective temperatures, and surface gravities of stars with given mass and chemical composition as a function of age. The tracks depend on the basic parameters such as the initial mass and chemical compositions. The properties of stellar populations in galaxies also depend on other properties not explicitly included in most of current stellar evolution models, i.e., the stellar rotation and close binary companions. The numerical calculations are also sensitive to the treatment of the convection, such as a mixing length parameter $`l/H_p`$ and convective overshooting.

The Padova stellar evolutionary tracks are adopted for stars from the main sequence (MS) to the early asymptotic giant branch (EAGB) stars. The Padova tracks cover wide ranges in age ($`0`$ to $`16`$ Gyr) and metallicity ($`Z=0.0001`$ to $`Z=0.05`$), and are one of the most complete sets of stellar evolutionary models currently available. The Padova tracks were calculated with revised radiative opacities (Iglesias, Rogers, & Wilson 1992) and with $`l/H_p=1.63`$. The adopted tracks were calculated for the following set of chemical compositions: ($`Y=0.230,Z=0.0001`$; Girardi et al. 1996), ($`Y=0.230,Z=0.0004`$; Fagotto et al. 1994a), ($`Y=0.240,Z=0.004`$; Fagotto et al. 1994b), ($`Y=0.250,Z=0.008`$; Fagotto et al. 1994b), ($`Y=0.280,Z=0.02`$; Bressan et al. 1993), and ($`Y=0.352,Z=0.05`$; Fagotto 1994a).

To obtain the isochrone for a given metallicity, the original isochrones are linearly interpolated in metallicity. Since the metallicity changes stellar evolution non-linearly, any extrapolation in metallicity should be avoided. Thus, stars of metallicity lower ($`Z<0.0001`$) and higher ($`Z>0.05`$) than the Padova tracks do not appear in our simulated CMDs and are not considered here. This is justified, because no significant populations of such extreme metallicities are known to exist in dwarf galaxies. In addition to these, close binary systems are not included, since the input tracks are only for isolated single stars. Although our code does not explicitly include binary stars, it effectively takes into account apparent binary stars as stellar blends. Gallart et al. (1999) reported a study that is more detailed on this issue.

For the RGB evolution, the mass loss law of Reimers (1977) is adopted;

$$\dot{M}=410^{13}\eta \frac{L/L_{}}{gR/R_{}}(M_{}\mathrm{yr}^1),$$
(8)
or equivalently,

$`\mathrm{log}(\dot{M})=`$ $``$ $`4.67+\mathrm{log}\eta +1.5\mathrm{log}(L/L_{})`$ (9)
$`\mathrm{log}(M/M_{})2\mathrm{l}\mathrm{o}\mathrm{g}T_{\mathrm{eff}},`$
where $`L,g,R,T_{\mathrm{eff}},`$ and $`M`$ are stellar luminosity, surface gravity, radius, effective temperature, and initial mass, respectively. The mass loss efficiency parameter $`\eta `$ is assumed to be $`1/3`$ (Renzini & Voli 1981).

The horizontal branch (HB) stars are distributed on the HR diagram according to a modified Gaussian mass distribution equation (Lee et al. 1990);

$`\mathrm{\Phi }(M)`$ $``$ $`\left[M(\overline{M_{\mathrm{HB}}}\mathrm{\Delta }M)\right](M_{\mathrm{RGB}}M)`$ (10)
$`\times \mathrm{exp}\left[{\displaystyle \frac{(\overline{M_{\mathrm{HB}}}M)^2}{\sigma ^2}}\right],`$
where $`\sigma `$ is a mass dispersion factor, and $`\overline{M_{\mathrm{HB}}}(M_{\mathrm{RGB}}\mathrm{\Delta }M)`$ is the mean mass of HB stars. A value of $`\sigma =0.06M_{}`$ has been chosen to represent the HB distributions similar to those of Galactic globular clusters (Rood 1973; Lee et al. 1990). This is equivalent to considering the statistical dispersions of mass loss rate along the RGB.

To convert the theoretical temperature-luminosity data to the observable colour-magnitude plane, a stellar spectral library by Lejuene, Cuisinier & Buser (1998) is used. This consists of Kurucz’s (1992) spectra for hotter stars (O$``$K), Bessell et al.’s (1989; 1991) and Fruks et al’s (1994) spectra for M giants, and Allard & Hauschudt’s (1995) for M dwarfs. In the original model spectra, systematic deviations become apparent when colour-temperature relations computed from the models are compared to the empirical ones at $`1Z_{}`$. The library adopted here is a version that the authors made by correcting the original library. The corrections are especially important for M star models. The fundamental parameters ($`T_{\mathrm{eff}}`$, $`g`$, and \[M/H\]$`=\mathrm{log}(Z/Z_{})`$) are wide enough to cover all spectral types and luminosity classes that appear in observed CMDs; $`3\times 10^3T_{\mathrm{eff}}5\times 10^4`$, $`1.02\mathrm{log}g5.5`$; $`5.0[\mathrm{M}/\mathrm{H}]1.0`$. To reduce computational time, synthetic magnitudes and colours tabulated by Lejuene et al. (1998) are linearly interpolated in log$`Z`$, log$`T_{\mathrm{eff}}`$, and log$`g`$.

### 3.2 Comparison with the CMD of the Draco dSph

In section 2, models B and C are shown to fit well the observed trends in the \[Mg/Fe\]$``$\[Fe/H\] diagram. Next, we discuss if the CMDs predicted by these models can reproduce the observed ones. Ideally, the SFH would be studied using a CMD which satisfies the following two conditions: (1) it covers a whole galaxy; (2) it is significantly deeper than the turnoff of the oldest stellar population. Unfortunately, such CMDs are not available for the galaxies discussed here. The situation, however, is better for the dSph Draco. Recently Piatek et al. (2001) imaged nine fields in and around the Draco dSph by using the KPNO 0.9m telescope and presented the CMDs down to a luminosity level $`2`$ mag fainter than the HB. Thus, we mainly discuss the dSph Draco in this section and briefly mention the dSphs Sextans and Ursa Minor.

Fig. 3 shows the observed CMD of the C0 field in the Draco dSph (Piatek et al. 2001), which is characterised by the narrow RGB and the well-populated red HB. The HB in figure 1 of Piatek et al. (2001) is somewhat distorted due to RR Lyrae stars which are marked with open circles in Fig. 3. Therefore, we do not take seriously the inconsistency between the observed and simulated blue HB morphology (V$``$R$`<0.1`$). Historically, Baade & Swope (1961) obtained the CMD of the Draco dSph and found that the RGB is generally similar to those of metal-poor globular clusters, although it is rather wide. The authors also found that the populous RHB is incompatible with the low metallicity of a Population II system. This is the notorious second-parameter problem of the dSph Draco.

Figs. 4a and 4b represent the CMDs simulated by models B and C, respectively, where the galactic age is assumed to be 12 Gyr, which is similar ages of Galactic globular clusters (e.g., Carretta et al. 2000); values of other parameters are written in table 1. The photometric errors and detection completeness in the simulations are taken from Piatek et al. (2001). Both models can reproduce characteristics of the CMD of the dSph Draco, such as narrow RGB and heavily populated RHB. Despite the metallicity dispersion ($`3[\mathrm{Fe}/\mathrm{H}]1.4`$), the theoretical CMDs results in narrow RGBs. This is due to the so-called age-metallicity degeneracy as explained below.

Fig. 5 shows CMDs simulated by assuming models B and C, where different colours represent different metallicity range as described in the figure. For presentation purpose, Figs. 5a and 5b contain larger numbers of stars than those of Figs. 4a and 4b, respectively. The adopted observational conditions are the same as those of Fig. 4. Clearly, there is little correlation between the colour and metallicity of RGB stars because of the age-metallicity degeneracy. Since stars of higher metallicity are younger, the metallicity and age effects on the RGB colours are cancelled out. This keeps the RGB narrow and tight. The age effect appears in the HB morphology too. Since only metal-poor and old ($`>10`$ Gyr) stars can evolve to BHB, younger core-He burning stars lie in the red part of the HB. The populous red HBs in Fig. 5 clearly show this age effect on the HB morphology. A weak correlation between the stellar colours and metallicities is found in the bright ($`m_V<18`$ or $`M_V<1`$) RGB.

In model C, the RHB is more populous and brighter than in model B. This is because the longer duration of the star formation period increases the number of younger and more massive core-He burning stars. However, the difference between models B and C is so small that it is difficult to determine the final epoch of star formation from the CMD that does not reach down to the turnoff level. Nevertheless, it can be said that, in the dSph Draco, chemical evolution was very slow and the duration of star formation was at least longer than 3.9$`6.5`$ Gyr.

This evolutionary picture can explain both the narrow RGB and the populous red HB of the Draco dSph. Thus, the second parameter problem in the dSph Draco is solved for the first time by introducing relatively long duration of star formation ($`>3.96.5`$ Gyr) which is fully consistent with the observed abundance patterns.

## 4 Discussion

The abundance patterns and the CMDs of the dSphs Draco, Sextans, and Ursa Minor share the common features. It can be said that they have similar SFHs and CEHs. The CMD of the Sextans dSph (see figure 2 in Suntzeff et al. 1993; figure 1 in Shetrone et al. 2001) is characterised by a narrow RGB and a red HB similar to the dSph Draco. This suggests that the SFH of the dSph Sextans is similar to that of the dSph Draco. In the CMD of the Ursa Minor dSph (see figure 1 in Shetorne et al. 2001), the blue HB is more populous and the RGB is narrower. This might suggest that the time scale of star formation was shorter in the dSph Ursa Minor. Nevertheless, it can be said that, in the Draco dSph, the SFR was very low and the duration of star formation was at least longer than 3.9$`6.5`$ Gyr. For the other dSphs, it is safe to say that the characteristics of the SFHs are the low SFRs and the long durations ($`>`$ several Gyrs) of the star formation period.

For the dSph Draco, model C may conflict with the CMD obtained by the Hubble Space Telescope (Grillmair et al. 1998). Through a comparison with the fiducial lines of metal-poor globular clusters M68 and M92, Grillmair et al. (1998) suggested that the dSph Draco is older than M68 and M92 by $`1.6\pm 2.5`$ Gyr. Their result, however, seems inconclusive. Firstly, the sparseness of the photometric sample makes it difficult to measure the contribution of any intermediate-age population. Although a predominant RHB is the characteristic of the dSph Draco, the few RHB stars appeared in their CMDs. Secondly, they used different techniques for photometry and presented the combined results. For bright stars aperture photometry was adopted, while the point spread function fitting method was used for faint stars. The border is unclear. This causes systematic discrepancies between the photometry of faint and bright stars. Third, a comparison with M68 and M92 is questionable. They used the F606W filter of the WFPC2 system. However, the WFPC2 Instrument Handbook recommends that the F555W filter is a better approximation to the Johnson V-band than the F606W. Their filter selection could have reduced the accuracy of the photometric calibration. Because of the combination of these uncertainties, it is dangerous to place much faith on their age estimation. More data are required to clarify these matters and images of a wide special coverage are particularly important.

In the above discussion, the decline in \[Mg/Fe\] observed in the dSphs is interpreted as the effect of the chemical enrichment by SNeIa. Since the signal-to-noise ratios of the data (Shetrone et al. 2001; S/N$`1336`$) for the dSphs are lower than other observations for metal-poor stars (e.g., McWilliams et al. 1995;S/N$`3040`$), one might argue that the quality of the data is not high enough to study the CEHs. However, the goal of Shetrone et al. (2001) was the measurement of overall abundance differences, not absolute values as presented in their study, and it can be said that \[$`\alpha `$/Fe\] declines at \[Fe/H\]$`2`$ on average (see figures 4 and 5 in Shetrone et al. 2001). Thus, it is a pertinent interpretation that enrichment by SNeIa causes the decline of \[$`\alpha `$/Fe\]. Spectroscopic observations of high signal-to-noise ratios are required to confirm this.

We have shown that low SFRs can explain the observed trends in the \[Mg/Fe\]$``$\[Fe/H\] diagram based on models assuming the Salpeter IMF. However, models of steeper IMFs and higher $`ϵ_{\mathrm{SF}}`$ can also reproduce the trends. The steeper IMFs lead to slow chemical evolution, so that SNe Ia start to explode at low metallicity. This results in the decline of \[Mg/Fe\].

Fig. 6 shows chemical evolution models E$``$G assuming different IMFs, where the IMF slopes are $`x=1.75,1.95,`$ and $`2.15`$, respectively, and $`ϵ_{\mathrm{SF}}=0.2`$ (a standard SFR for the solar neighbourhood; Arimoto et al. 1992). These models, except for model E, are also consistent with the observations. As with model A, the abundance ratios predicted by model E are too high and thus inconsistent with the observations. For the observed metallicity range (\[Fe/H\]$`<1.4`$), the trends of \[Mg/Fe\] of models E$``$G are very similar to those of models A$``$D, although they become different at higher metallicity (see Fig. 1).

Figs. 7a and 7b represent the CMDs simulated by models F and G, respectively. The photometric errors and detection completeness are taken from Piatek et al. (2001), i.e., the same as the simulations of models B and C shown in Fig. 4. The CMDs in Fig. 7 resemble those in figure 4 in terms of the morphologies of RGB and HB. The slow chemical evolution and the long ($`>1.6`$ Gyr) durations of star formation of models F and G result in a narrow RGB and the populous RHB by the same reasons as those of models B and C. This indicates that both abundance patterns and morphologies of the RGB and HB become similar when models have similar age-metallicity relations (AMRs).

Figs. 8a and 8b show the AMRs of models A$``$D and E$``$G, respectively. The figures demonstrate that the AMRs of models E$``$G with the steeper IMFs result in very similar to those of models A$``$D with the lower SFRs at low-metallicity (\[Fe/H\]$`<1.4`$). Since both the abundance patterns and the morphologies of the RGB and HB are affected by CEHs, an additional information is necessary to solve the degeneracy of IMF and SFR. In the models E$``$G the initial star formation lasts only $`1.62.2`$ Gyrs, while in the models B$``$C it continues as long as $`3.96.5`$ Gyrs. Therefore, the models B$``$C should give systematically brighter turnoff magnitudes than the models E$``$G, and magnitude difference between the main sequence turnoff and the HB $`V(\mathrm{TO}\mathrm{HB})`$ in the models E$``$G should be smaller than in the models B$``$C. The models E$``$G give $`V(\mathrm{TO}\mathrm{HB})`$ is $`3.1`$ mag, while the models B and C does $`V(\mathrm{TO}\mathrm{HB})3.1`$ mag and $`2.8`$ mag, respectively. Deeper images with wider field of views will allow us to determine the final epochs of the initial star formation.

Figs. 1 and 6 show that \[Mg/Fe\] converges at higher metallicity in the models with the different SFRs and that it does not in those with the different IMFs. This difference allows us to derive the IMF and SFR independently. Abundance patterns in more massive dwarf galaxies and/or dwarf irregular galaxies will test whether the IMF or SFR produces in the abundance patterns and low metallicities in the dSphs. We note that even steeper ($`x>2.75`$) IMFs are required to explain the low metallicities of dwarf irregular galaxies if the IMF is the primary cause of the low metallicities. The CMD of the dwarf irregular galaxies in the Local Group reveal that they contain very old ($`>10`$ Gyr) populations as well as young ones, i.e., the star formation and chemical enrichment have continued for at least 10 Gyr. Since they are still metal-poor, extremely steep ($`x>2.75`$) IMFs are required, which have never been reported. Therefore, we believe that low SFRs are the primary characteristic of the SFHs of the dSphs.

So far, the gas infall and outflow during star formation have been neglected, since no evidence of the gas infall and outflow has presented. Metallicity distribution could provide a clue to judge whether gas infall or outflow should be considered. We stress that spectroscopic observations is indispensable to derive a metallicity distribution. Fig. 5 demonstrates little correlation between the colours and metallicities of the RGB stars. Because of the age-metallicity degeneracy and contamination by AGB stars, a stellar colour is not always a good indicator of the metallicity. Thus, a careful procedure should be adopted to convert stellar colours into metallicities (Ikuta 2001), although a relatively simple technique is often used (e.g., Harris & Harris 1999).

Long lasting star formation at a very low rate explains the observed trend of \[Mg/Fe\] and the CMDs. Because of the low SFRs, plenty of gas ($`97`$ percent of the galaxy mass) still remain even at the final epoch of star formation (i.e., $`t=\mathrm{\Delta }T_{\mathrm{SF}}`$). The remaining gas has to be removed to complete star formation and to evolve to a gas-poor system. Energy injection from supernovae could not be a sufficient mechanism to expel the gas from dwarf galaxies (MacLow & Ferrara 1999; Ikuta 2001). If this is the case, gas removal from dwarf galaxies should result from external mechanisms such as ram pressure stripping and/or tidal shocks. Studying the SFHs of 18 dwarf galaxies in the Local Group based on the CMDs derived by a uniform method of photometry, Ikuta (2001) found that the durations of star formation period correlate with the distances from the Milky Way or M31. The correlation may suggest that the environmental effects play a key role in the evolution of the dwarf galaxies. Numerical simulations (Moore et al. 1998; Mayer et al. 2001) clearly showed that late-type dwarf galaxies entering the dark matter halo of a massive galaxy are transformed into early-type owing to repeated tidal stripping and dynamical instabilities. Thus, we conclude that the star formation in the dSphs Draco, Sextans, and Ursa Minor terminated due to gas stripping by the Milky Way.

The hierarchical clustering galaxy formation model suggested that star formation in dwarf galaxies can only occur before the cosmological re-ionisation epoch and there is no major star formation activities later because gas cannot cool and condense to form stars (e.g., Cen 2001). This appears to be inconsistent with stellar populations observed in the Local Group dwarf galaxies. CMDs obtained by the HST (e.g., Ikuta 2001) and our study presented in this paper clearly show extended and recent star formation activities in the Local Group dwarf galaxies.The issue is still being debated. However, a more recent simulation (Kitayama et al. 2001) suggested that star formation can occur in small objects if they have baryonic mass larger than the threshold mass of $`10^9M_{}`$ and $`10^6M_{}`$ at redshifts of $`<3`$ and $`5`$, respectively. The masses of today’s dSphs in the Local Group are $`10^710^8M_{}`$ (e.g., Mateo 1998) and our results imply that progenitors of today’s dSphs lost $`97`$ % of their mass after the long lasting star formation. Therefore, the masses of the progenitors are estimated at $`10^910^{10}M_{}`$. Since these exceed the threshold predicted by Kitayama et al. (2001), the SFHs derived here do not contradict recent simulations of the formation of dwarf galaxies.

## 5 Conclusion

Based on comparisons between the theoretical chemical evolution models and the observed abundance patterns, we conclude that the initial star formation continued for a long duration ($`>3.96.5`$ Gyr) in the dSphs Draco, Sextans, and Ursa Minor. Our simulation of the CMDs shows that the long duration of star formation can solve the second parameter problem of the Draco dSph. Because of the age-metallicity degeneracy, the RGBs are kept narrow and tight despite their large metallicity dispersions.

We have discussed the two cases which are consistent with the observed abundance patterns and the CMDs. The first case is a combination of low SFRs ($`ϵ_{\mathrm{SF}}=110^2510^3`$ Gyr-1) and the Salpeter IMF ($`x=1.35`$), while the second is a combination of the solar neighbourhood SFR ($`ϵ_{\mathrm{SF}}=0.2`$ Gyr-1) and the steeper IMFs ($`x=1.752.15`$). The two cases are discriminated neither by the abundance patterns nor by the CMDs of bright stars. However, the initial star formation period is as long as $`3.96.5`$ Gyrs if the SFR is low, while it is much shorter ($`1.62.2`$ Gyrs) if the SFR is high. Thus, turnoff magnitude should be different between the two cases. Deeper images with wider field of views allow us to determine the final epochs of the initial star formation.

###### Acknowledgements.

We wish to express my gratitude to the anonymous referee for very helpful suggestions and comments. We are also grateful to T. Kodama, and H. Susa, for fruitful discussion and comments and to B. Jones who carefully read this manuscript and gave us helpful comments and suggestions. C.I. wishes to thank the Japan Society for Promotion of Science for financial support.


text:
# 1 The KNO multiplicity distribution, 𝑃(𝑛)=𝜎_𝑛/𝜎_{𝑡𝑜𝑡}. The points are the E735 Tevatron data []. The range A corresponds to the multiperipheral kinematics. The range B corresponds to the non-interacting, ideal, gas approximation, 𝑛≪𝑛_{𝑚𝑎𝑥}. C is the VHM region.

Phenomenology of Very High Multiplicity Hadron Processes
J.Manjavidze111JINR,Dubna, Russia $`\&`$ Inst. of Phys., Tbilisi, Georgia. E-mail: joseph@nusun.jinr.ru and A.Sissakian222JINR, Dubna, Russia. E-mail: sisakian@jinr.ru

## Abstract

We discuss the possibility to suppress the nonperturbative effects if the very high multiplicity hadron final states are chosen. The theoretical uncertainties and possible experimental measurements are described.

1. Very high multiplicity (VHM) events are expected to be produced at the future experiments at BNL (STAR), Fermilab (CDF), CERN (ATLAS). Then, it is natural to expect the VHM region to give new important information, which will be complementary to the ongoing investigations and/or unavailable in other than the VHM process. The studies of the VHM region have been proposed recently in  and attract an increasing interest in high energy physics .

In this Letter, where we briefly consider the VHM topic in whole, we particularly emphasize the related physical quantities expected to be measured with forthcoming experiments. In this sense, it is crucial to understand experimental abilities of such measurements, since the topological cross section, $`\sigma _n`$, is extremely small (e.g., at TeV energies, $`\sigma _n10^7\sigma _{tot}`$, compared to the total cross section $`\sigma _{tot}`$, Fig.1) and could lead to certain problems.

Started investigations of VHM processes, one has first of all ask a question: why the mean multiplicity in hadron processes being energy dependent as $`\mathrm{ln}^2(s/m_\pi ^2)`$ is much $`smaller`$ than the threshold multiplicity value $`n_{max}=\sqrt{s}/m_\pi `$ ? Here $`\sqrt{s}`$ is the c.m.s. energy and $`m_\pi `$ is the pion mass.

In the present Letter, we follow the idea  that the multiple production process could be considered as a process of hot tea cooling in a cold environment. Then very hot tea cooling in very cold room is analogous to the VHM process. In that case, the cooling process proceeds "quickly", i.e. any process obstacled to the cooling fluctuations must be suppressed. The idea that the hadron multiple production process may consist of various mechanisms has been considered in .

In the QCD parton picture  the "fastness" of a process means that parton decay must prevail over its dissipation from an interaction zone. This is possible if the parton virtuality $`|q|`$ is high enough as far as the parton lifetime is $`1/|q|`$. As a result, the VHM process should be hard and, therefore, the forces preventing particle production (confinement) must be negligible.

Fig.2 shows the longitudinal-transverse momenta phase space domain of the VHM process. It is remarkable that the VHM domain occupies both regions, the "low-$`p_{}`$" and the "low-$`p_{}`$" ones. Therefore, despite the fact that the VHM process is hard one can not use the Leading Logarithm Approximation (LLA)  to describe it . To stress is that, following the LLA formalism, there are the only two domains, namely the Regge and Deep inelastic Scattering (DIS) domains, as shown in Fig.2. The corresponding cross sections may be sufficiently large and the LLA is applicable.

The "low-$`p_{}`$", or "low-$`x`$"  in DIS, also means that VHM particles have low energies. For that reason, the hadrons are mainly resulted from the "gluing" of coloured constituents created in the hard process. Therefore, the (soft) process of coloured parton pair production from the vacuum is negligible in the VHM domain. This distinguishes VHM processes from other hard processes (e.g., from hadron production in $`e^+e^{}`$-annihilations, or DIS processes). The VHM dynamics is free of the non-perturbative background.

Sec.2 argues the above conclusion, in Sec.3 the main properties of VHM dynamics are discussed and the concluding remarks are given in Sec.4.

2. The dominance of "Regge" processes stands for the "softness" of hadron dynamics . This is based on the experimental observations: the mean transverse momentum of secondary particles, being energy and multiplicity independent, has a limited value, and the diffraction radius of hadron elastic scattering increases with increasing energy. Such "multiperipheral" kinematics is described by the Regge model with the soft Pomeron ansatz.

From the "microscopical" point of view, it means that the soft process can be roughly described by the very spatial class of the Feynman "ladder" diagrams of perturbative QCD (pQCD). This is the so-called "BFKL Pomeron" solution . It is based on the LLA approach and assumes that the process described is sufficiently hard333The point is that this Regge-like theory depends on the phenomenological energy scale parameter $`s_0`$ and, therefore, its range of application may be chosen arbitrary..

The both models are based on the assumption that the mean transverse momentum of produced hadrons is small. In the frame of the soft Pomeron approach, this constrain is hidden in the conservation laws, which are the consequence of the underlying non-Abelian gauge symmetry. In the BFKL Pomeron case, the LLA place a role of the constrain.

It is worse to note that the constrain restricts the production dynamics and, therefore, the produced states can not entirely cover the phase space with the same density, Fig.2. This requires the constrain to be suppressed in order to lead to the VHM final states. Remembering that the constrains are the consequence of the long-range connections among interacting degrees of freedom, this also points out that the VHM process must be the hard process.

The VHM scenario discussed realizes if (i) the hard processes has place in the VHM region and if (ii) the hard process ensures the "fast" hadron production.

The formal proof of the item (i) contains the following steps . First, it can be shown that if the final-state interactions are excluded, then the $`soft`$ process only satisfies the following asymptotics:

$$\sigma _n^s<O(e^n),$$

i.e. the "soft" topological cross section $`\sigma _n^s`$ falls down faster than any power of $`\mathrm{exp}\{n\}`$. The reason is just the softness of the process. On the other hand, the $`hard`$ process leads to the following asymptotic estimate:

$$\sigma _n^h=O(e^n).$$

Consequently, one always finds such energy $`\sqrt{s}`$ and multiplicity $`n`$ that

$$\sigma _n^h\sigma _n^s\mathrm{for}n\sqrt{s}/m_\pi .$$

This proves the item (i).

The qualitative argument for the statement (ii) follows from the Ehrenfest-Kac model  of the irreversibility phenomenon444The Ehrenfest-Kac model considers two box $`a`$ and $`b`$ and $`2N`$ $`numbered`$ balls. One may assume that the ”far from the equilibrium” initial state presents the situation when all $`2N`$ balls are e.g., in the box $`a`$ and in the equilibrium each box contains $`N`$ balls. Then, choosing $`randomly`$ the number of balls, one should take the corresponding ball from the box and put it into the other one. At the following step, the number of a ball will be again randomly chosen, and the corresponding ball should be put into the other box, and so on. One may measure the time $`t`$ in the units of the ball shifting time. Then, the number of the ”produced particles” in the box $`b`$ measures the trend to equilibrium. In this equilibrium state, $`N_a=N_b`$ (Fig3).. It shows that if the initial state is far from the final one, then the system undergoes to equilibrium in the fastest way 555This phenomenon was not mentioned in the known publications concerning the Ehrenfest-Kac model.. The result of the computer modelling of "particle production" within discussed model, Fig.3, confirms this. Indeed, as it is seen there are no fluctuations at early stage of the process: number of the "produced particles" is proportional to $`t`$, up to $`t1000`$.

The formal proof of the statement (ii) can be obtained from the KNO-scaling of the multiplicity distribution in pQCD jets . Introducing the quantity

$$\mu (n,s)=\frac{1}{n}\mathrm{ln}\{\frac{\sigma _n(s)}{\sigma _{tot}(s)},\}$$
(1)
one finds at the mass $`M`$ of a jet,

$$\mu _1(n,M)=\frac{a}{\overline{n}_j(M^2)}+O(1/n).$$

Here $`a>0`$ and $`\overline{n}_j(M^2)`$ is the mean hadron multiplicity in a QCD jet, $`\mathrm{ln}\overline{n}_j(M^2)\{\mathrm{ln}(M^2/\mathrm{\Lambda }^2)\}^{1/2}`$. The two-jet contribution at the same energy $`M`$ reads as

$$\mu _2(n,M)=\frac{a}{\overline{n}_j(s/4)}+O(1/n)>\mu _1(n,s),$$

where the energy conservation law is taken into account. Therefore, the asymptotics over $`n`$ is defined by the events with the lowest numbers of pQCD jets.

This confirms the statement (ii) since just the jet mechanism of particle production is known to be the fastest one. Moreover, we conclude that the number of the jets decreases with increasing total multiplicity $`n`$ and, thus, the heavier jets are produced in the VHM region. However, one has to note that this conclusion is based on the LLA and, therefore, we consider it just as a qualitative prediction.

3. From the above, one concludes that the VHM process is a hard process. In this case, the process of incident energy dissipation evolves $`freely`$, i.e. without the influence of the "hidden constrains". This is the important observation because it allows to assume that in such conditions the of final state distribution trends to the "equilibrium".

Let us remind some statistical properties of the maximally constrained systems, also called the "exactly integrable systems". It is known that if the system has $`2N`$ degrees of freedom and holds $`N`$ constrains (first integrals in involution) then no thermalization occurs in such a system. The first observation of this phenomenon belongs to Fermi, Pasta and Ulam . On the other hand, it is evident that the system without any constrains should drift to the equilibrium.

The hadron dynamics hides definite number of constrains, which is not enough to completely suppress particle production, i.e. the thermalization process. This conclusion follows from the existence of large multiplicity fluctuations, see Fig.1. Nevertheless, the constrains occur and the mean multiplicity is reasonably small $`\overline{n}(s)n_{max}`$.

To stress is that the term "equilibrium" is understood here as the absence of any macroscopic correlations among the extensive parameters of a system. For example, the absence of energy correlations implies the thermal equilibrium. Then, one can obtain the experimental condition, which reflects the "equilibrium phenomena" . It has been shown that the condition

$$R_l(n,s)=\frac{|K_l(n,s)|^{2/l}}{|K_2(n,s)|}<1,l=3,4,\mathrm{},$$
(2)
is the necessary and sufficient one to reach thermalization. Here $`K_l(n,s)`$ are the $`l`$-point energy correlators in the $`n`$-particle system, $`l<n`$. To note is that (2) reminds the "correlation relaxation condition" of Bogolyubov .

The tendency of the system to reach the equilibrium is considered as the dynamical effect, when the (symmetry) constrains are switched off, and for that reason the system may drift to the equilibrium state.

Fig.4 shows the result of Monte Carlo PYTHIA event generator  simulation of the ratio $`R(n,s)`$ for the domain A in Fig.1. The absence of the thermalization process is natural if one notes that PYTHIA resembles the Regge model. At the same time, it can be shown that in the region B the system without fail reaches the thermal equilibrium, $`R(n)1/n`$ in that domain.

4. We can conclude that:

(i) the VHM process allows us to investigate the hadron dynamics beyond the standard multiperipheral kinematics;

(ii) it is the hard process in sense that the influence of hidden constrains, connected with the high symmetries and the nontrivial vacuum of the Yang-Mills theory, is not important in the particle production dynamics;

(iii) the VHM process allows us to reach the thermal equilibrium, where the ordinary thermodynamical methods can be used. (This extends experimental possibilities and, therefore, is extremely important. For example, the $`\mu `$-quantity introduced in (1) is the "chemical potential" if measured in the units of the mean energy of produced particles );

(iv) the VHM domain seems to be one of the mostly interesting regions to investigate the collective phenomena, e.g. phase transitions. (To note is that the VHM final state is a "cold" one).

To note is that the production dynamics of VHM states might be well enough described by pQCD, but the logarithmic accuracy seems to be not enough. This is a main theoretical problem and to this end a more accurate perturbation theory, described in , has been formulated.

At the end, let us add that, as it follows from the above conclusion, the VHM process leads to the $`dense`$, $`cold`$ and $`equilibrium`$ locally-coloured state, e.g. "cold" plasma 666It is hard to imagine the coloured state to form the QCD plasma one if the former is ”hot”. In this case the kinetic motion would rapidly separate the charges and this must lead to the strong polarization of the vacuum. This effect is absent in the QED plasma..

Acknowledgements

Authors would like to take the opportunity to thank Yu. Budagov, V.Kadyshevski, A.Korytov, E.Kuraev, L.Lipatov, V.Matveev, V.Nikitin and E.Sarkisyan for fruitful discussions and constructive comments. Special thanks to STAR (BNL) and CDF (FNAL) physics communities for continuous interest in the VHM problem.