format: text

text:
# Interplay of quantum and thermal fluctuations in a frustrated magnet

## I Introduction

In geometrically frustrated magnets, the arrangement of the spins on a lattice precludes the establishment of a simple collinear Néel state. In the absence of such an ‘obvious’ ordering pattern, these magnets are often subject to a large number of competing instabilities. This endows them with a rich behaviour, manifested in unconventional phases, phase transitions and excitations, a large low-energy density of states and unusual dynamical properties.reviews

When added to a frustrated magnet which is not ordered classically at temperature $`T=0`$, both thermal and quantum fluctuations often generate an ordering transition, a process known as order by disorder.villainobdo ; shenderobdo ; henleyobdo  Whereas quantum order by disorder occurs in a fashion quite analogous to the ordering induced by thermal fluctuations, the nature of the ordering (if any) of even some ‘simple’ model systems (such as the Heisenberg magnet on the pyrochlore lattice) in the presence of stronger quantum fluctuations remains unclear.

The interplay of thermal and quantum fluctuations has so far not received a great deal of theoretical attention for frustrated magnets. By contrast, there has been a considerable amount of interest devoted to the real-time dynamics and transport at finite temperature near quantum critical points.kedar  In addition, this combination has even been claimed to provide an alternative use of quantum effects for speeding up computations.appros

In this publication, we study the combined effect of thermal and quantum fluctuations on possibly the simplest quantum frustrated magnet, namely an Ising model on the triangular lattice in a transverse magnetic field. The manifold of classical ground states is hugely degenerate, and correlations averaged over it are critical.wannhout  This model thus enables us to investigate how the two types of fluctuations together establish and destroy order out of the exponentially large set of classical ground states. By changing the values of transverse field, $`\mathrm{\Gamma }`$, and temperature, $`T`$, their strengths are in principle separately experimentally tunable.

We map out the phase diagram of the Ising model on the triangular lattice in the transverse field-temperature plane. The bulk of the paper is devoted to a numerical study based on a continuous time quantum Monte Carlo algorithm. We find that there are three different phases. In the regime where quantum fluctuations are weak, yet dominate over thermal fluctuations, they generate a three-sublattice ordered phase. Upon increasing the strength of thermal fluctuations, it ‘melts’ into a critical phase which has a finite extent and displays drifting exponents. This dome-shaped critical phase is delineated above and below by Kosterlitz-Thouless (KT) transitions at finite temperature. At $`T=0`$, it terminates in the classical critical point at zero field and in a quantum phase transition at a finite $`\mathrm{\Gamma }_c`$. At high temperatures and large field strengths, one retrieves a simple paramagnet – order by disorder is a reentrant phenomenon. We briefly comment on the advantages of different diagnostics for the floating KT phase.

The results reported here largely build on the following pieces of work. In Ref. ftfim, , a number of frustrated transverse field Ising models were studied, and a connection was pointed out between these models and classical stacked magnets in a scaling limit. An influential paper by Blankschtein and coworkersblankschtein  on stacked triangular magnets derived a Landau-Ginzburg-Wilson theory for these systems, and predicted the stacked problem to be in the universality class of the $`XY`$ model in three dimensions with a six-state clock term, a problem which has in turn been studied by Jose et al.Jose77  Collecting these results together led to a conjecture of the phase diagram we map out in the following.ftfim  We note that, by adding a continuous degree of freedom to a classical triangular Ising antiferromagnet, a floating critical phase has previously been obtained by Chandra et al.premi94

The general validity of the Landau-Ginzburg approach has been somewhat called into question by Monte Carlo simulations on the ferromagnetically stacked triangular Ising antiferromagnet, with several groups presenting evidence in favour of critical behaviour in the three-dimensional universality class, while others found a different, possibly new, universality class to be a more likely scenario.heino89 ; tricrit ; bunker ; plumer95 ; collinsrev  Our work here lends further support to the Landau-Ginzburg approach which is independent of the previous numerical simulations.

## II Model

Possibly the simplest model which has both a high classical ground state degeneracy and a non-trivial quantum dynamics is the transverse magnetic field antiferromagnetic Ising model on the triangular lattice. The model has the following Hamiltonian:

$$H=J\underset{i,j}{}\sigma _i^z\sigma _j^z\mathrm{\Gamma }\underset{i}{}\sigma _i^x,$$
(1)
where $`J>0`$ is nearest-neighbour exchange couplings and $`\mathrm{\Gamma }`$ is the magnetic field strength.

The classical model triangular Ising magnet is obtained in the absence of a transverse field, $`\mathrm{\Gamma }=0`$. Its solution is well known:wannhout  it is disordered at any finite temperature, and critical at $`T=0`$, where it retains an extensive zero point entropy, $`𝒮`$, with $`𝒮/N=0.323k_B`$, where $`N`$ is the number of spins and $`k_B`$ denotes the Boltzmann constant. Upon application of an infinitesimal transverse field, $`\mathrm{\Gamma }=0^+`$, at $`T=0`$, the magnet orders into a three sublattice $`\sqrt{3}\times \sqrt{3}`$ pattern with a sublattice magnetisation pattern $`(1,0,1)`$.ftfim

### II.1 Mapping to a stacked magnet

Using the Suzuki-Trotter formalism,Suzuki71 ; Trotter59  one can map the d-dimensional quantum model onto the (d+1)-dimensional classical problem. The partition function of the quantum Hamiltonian (Eq. 1), $`\mathrm{Tr}\mathrm{exp}(\beta H)`$, is equivalent to the following partition function:

$`Z=\mathrm{Tr}\mathrm{exp}(S_{cl}),`$
$`S_{cl}=K{\displaystyle \underset{i,j,\tau }{}}S_{i,\tau }S_{j,\tau }K^\tau {\displaystyle \underset{i,\tau }{}}S_{i,\tau }S_{i,\tau +1},`$ (2)
where $`S_{i,\tau }`$ are classical Ising spins taking values $`1`$ or $`1`$; $`\tau `$ is the index in the imaginary time direction and runs from $`0`$ to $`N`$, $`N=\beta /ϵ`$ is the number of slices in the imaginary time direction and $`ϵ`$ is the imaginary time step; $`K=Jϵ`$ and $`K^\tau =\frac{1}{2}\mathrm{ln}\mathrm{tanh}ϵ\mathrm{\Gamma }`$. The mapping becomes exact in the scaling limit $`ϵ0,K^\tau \mathrm{}`$, while maintaining

$$e^{2K^\tau }=ϵ\mathrm{\Gamma }.$$
(3)
It is in this formulation that the dimensional crossover between the high temperature $`d=2`$ and the low temperature $`d=3`$ behaviour is most transparent. Besides providing an efficient Monte Carlo algorithm (see below), the above mapping also allows us to make contact with the literature on stacked triangular magnets,collinsrev  from which our problem can be obtained in the scaling analysis of infinite coupling strength in the third dimension, with the third dimension being of finite length at nonzero temperature.

### II.2 LGW analysis

Blankschtein et al.blankschtein  have identified a complex XY order parameter for a stacked triangular antiferromagnet: the amplitude and phase of the Fourier mode at $`\stackrel{}{Q}_\pm =(\pm 4\pi /3,0)`$, $`\psi _\pm =m\mathrm{exp}(\pm i\varphi )`$. They have obtained the following LGW Hamiltonian (see below):

$`H_{LGW}={\displaystyle \underset{\stackrel{}{q}}{}}(r+q^2)m^2+u_4m^4`$
$`+u_6m^6+v_6m^6\mathrm{cos}(6\theta ).`$ (4)
A six-fold clock anisotropy thus appears at sixth order in $`m`$. The six-fold clock term is irrelevant at $`d=3`$ so that the transition into a phase with nonzero $`m`$ should be in the $`d=3`$ XY universality class. This applies both to the zero temperature quantum phase transition as we vary $`\mathrm{\Gamma }`$ and to the corresponding transition in a stacked Ising magnet. With the clock term being dangerously irrelevant, the transition is immediately into a state with a broken clock symmetry, the details of which depend on the sign of $`v_6`$ (see below).

At finite temperatures the triangular transverse field Ising model maps onto a stacked magnet of finite size, $`L_\tau `$, in the temporal direction, which is effectively $`2`$-dimensional as the temporal correlation length, $`\xi _\tau `$, exceeds $`L_\tau `$ close to a continuous phase transition. The $`2`$-dimensional six-state clock model, studied by Jose et al.,Jose77  has a remarkable phase diagram consisting of three phases: a disordered phase, an extended intermediate KT phase, and an ordered phase. The extended KT phase owes is existence to the fact that in $`d=2`$, the clock term does not become relevant immediately below the KT transition of the XY model; instead, the coupling needs to be increased further until the clock term is sufficiently strong to assert itself.

One of the central objectives of this paper is to establish the presence of this floating KT phase as a finite-temperature induced dimensional crossover in the $`2+1`$ dimensional quantum Ising magnet on the triangular lattice. In our context, the presence of such a phase would be the result of an XY order parameter generated by frustration and stabilised by quantum fluctuations, together with destabilising fluctuations provided by the finite temperature.

Details of the ordered state depend on the sign of $`v_6`$. $`H_{LGW}`$ is minimised by $`M=|M|\mathrm{exp}(i\mathrm{\Phi })`$, with $`\mathrm{\Phi }=n\pi /3`$ and $`\mathrm{\Phi }=(n+1/2)\pi /3`$ for $`v_6<0`$ and $`v_6>0`$ respectively, with $`n=1\mathrm{}6`$. The real space configurations of the ordered phases can be obtained by Fourier transforming these modes. In the case of $`v_6<0`$, Fourier transforming yields the following magnetisations assigned to the three sublattices of the triangular lattice: $`|M|(1,1/2,1/2)`$. We depict this phase as $`(+)`$. In the case of $`v_6>0`$, we have the following sublattice magnetisations: $`|M|(\sqrt{3}/2,0,\sqrt{3}/2)`$. We depict this phase as $`(+0)`$. There are six degenerate states for both $`v_6<0`$ and $`v_6>0`$.

It is difficult to determine the sign of $`v_6`$ – in particular, in the effective Hamiltonian, its value may drift and even change sign.blankschtein ; miyashita  The partition function $`Z=\mathrm{exp}(_{ij}s_iK_{ij}s_j)`$ of the Ising model $`H=(1/\beta )_{ij}s_iK_{ij}s_j`$ is equivalent to the following partition function (see, e.g. Ref. amit, )

$`Z{\displaystyle }D\psi \mathrm{exp}()={\displaystyle }D\psi \mathrm{exp}\{{\displaystyle \underset{ij}{}}\psi _iK_{ij}\psi _j`$
$`+{\displaystyle \underset{i}{}}\mathrm{ln}\mathrm{cosh}\left({\displaystyle \underset{j}{}}2K_{ij}\psi _j\right)\},`$ (5)
The first term is the interaction matrix, which determines the soft modes $`\psi _\pm `$ to be located at $`Q_\pm `$ (Eq. 7). Performing the Fourier transform and expanding $`\mathrm{ln}\mathrm{cosh}(x)=\frac{x^2}{2}\frac{x^4}{12}+\frac{x^6}{45}\frac{17x^8}{2520}+\mathrm{}`$, we obtain the following expression for the contribution of the non-linear term to $``$

$`_{nl}`$ $`=`$ $`2{\displaystyle \underset{𝐤}{}}K(𝐤)K(𝐤)\psi (𝐤)\psi (𝐤)`$
$`+`$ $`{\displaystyle \frac{4}{3}}{\displaystyle \underset{𝐤_i}{}}\left[{\displaystyle \underset{i=1}{\overset{4}{}}}\left(K(𝐤_i)\psi (𝐤_i)\right)\delta ^{}(𝐤_1+\mathrm{}+𝐤_4)\right]`$
$``$ $`{\displaystyle \frac{64}{45}}{\displaystyle \underset{𝐤_i}{}}\left[{\displaystyle \underset{i=1}{\overset{6}{}}}\left(K(𝐤_i)\psi (𝐤_i)\right)\delta ^{}(𝐤_1+\mathrm{}+𝐤_6)\right]\mathrm{}.`$
The interactions for the stacked triangular antiferromagnet are written in Fourier representation as

$`K(𝐤)={\displaystyle \underset{k_x,k_y,k_z}{}}(J[\mathrm{cos}(k_x)+\mathrm{cos}(k_x/2+\sqrt{3}k_y/2)`$
$`+\mathrm{cos}(k_x/2\sqrt{3}k_y/2)]J^{}\mathrm{cos}(k_z)),`$ (7)
where the sum is over the hexagonal Brillouin zone. The LGW Hamiltonian is thus constructed in terms of the coefficients $`m\mathrm{exp}(\pm i\varphi )`$), varying slowly in space, of the ‘soft’ modes $`\psi _\pm `$.

The primed delta functions indicate that the wavevectors $`𝐤_i`$ need to add up to a reciprocal lattice vector. This requirement leads to the XY nature of the effective theory at low order. At sixth order, $`\psi _+^3\psi _{}^3`$ and $`\psi _+^6+\psi _{}^6`$ occur, the latter being the clock term, $`m\mathrm{cos}(6\theta )`$. $`v_6<0`$ at this order. However, the $`\mathrm{cos}(6\theta )`$ term appears at higher orders as well, e.g., $`\psi _+\psi _{}(\psi _+^6+\psi _{}^6)`$ at eighth order. Since the terms in the series have alternating sign and large coefficients, one cannot reliably determine the sign of $`v_6`$ near the transition this way.

## III Monte Carlo method

### III.1 The quantum Monte Carlo algorithm

The $`d=2+1`$ dimensional classical problem obtained this way has only positive weights – there is no ‘sign problem’. It can therefore be studied reasonably straightforwardly by Monte Carlo simulations. We thus simulate the classical problem defined by Eq. 2. However, it is difficult to simulate the discretized version of Eq. 2 because of the scaling limit (Eq. 3): in order to avoid discretisation errors, one has to take a very large ferromagnetic coupling in the imaginary time direction $`K^\tau \mathrm{}`$ at the same time increasing the height of the system exponentially as $`e^{2K^\tau }`$.

In order to avoid this problem, we use a continuous time algorithm.beardwiese ; prokofev ; Rieger99  The basic idea behind this algorithm is that, in the scaling limit, the density of domain walls in the imaginary time direction becomes exponentially sparse and it is thus more efficient to keep track of the location of the domain walls, using $`\mathrm{exp}(2K^\tau )`$ as a unit of length. Thus, the height of the system, $`N/\mathrm{exp}(2K^\tau )=\beta /(ϵ\mathrm{exp}(2K^\tau ))=\beta \mathrm{\Gamma }`$, measured in units of $`\mathrm{exp}(2K^\tau )`$, remains fixed in the continuum limit.

Due to the frustrated nature of our problem we cannot use a cluster algorithmSwendsen87 ; Wolff89 ; Krauthreview  in the space directions as one could for the case of unfrustrated magnets; however, we can use a cluster algorithm in the time direction. The algorithm works as follows (for more details, see Ref. Rieger99, ). We pick a random site on the triangular $`L\times L`$ lattice and build a cluster on that site in the imaginary time direction. Its length, $`\tau `$, is distributed according to the probability distribution $`P(\tau )\mathrm{exp}(\mathrm{\Gamma }\tau )`$ – this prescription eliminates a ‘freezing’ of the algorithm due to the diverging temporal coupling. A given cluster is flipped using the Metropolis prescription in the spatial direction, i.e. we flip the cluster with probability $`p=\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }E))`$, where $`\mathrm{\Delta }E`$ is the (spatial) energy difference between the original configuration and the configuration with a flipped cluster.

One check of our Monte Carlo algorithm consisted of comparing it with the diagonalisation of a $`3\times 3`$ lattice, and we have found excellent agreement between the two approaches.

### III.2 Parameters of the simulations

We impose periodic boundary condition on the triangular lattice and performed simulations on lattices of size $`L=9,12,15,18,24,30,36,48,60,75`$. We have estimated the correlation time $`\tau _c`$ for different parameters. We perform usually $`500\tau _c`$ Monte Carlo steps (MCS) for equilibration and from $`10^4\tau _c`$ to $`510^5\tau _c`$ MCS for averaging.

## IV Scaling analysis

In this section, we outline the quantities useful for detecting a possible floating KT phase and the adjacent KT phase transitions as well as the ordered phase.

### IV.1 Order parameters

The complex XY order parameter deduced from the LGW analysis is given by:

$$me^{i\theta }(m_1+m_2e^{i\frac{4\pi }{3}}+m_3e^{i\frac{4\pi }{3}})/\sqrt{3}$$
(8)
where the $`m_i`$ are the magnetisations of the three sublattices, and $`m`$ is real and positive. $`m`$ is close to 1 in the limit of zero temperature and vanishing transverse field. It equals zero for the disordered and KT phases, vanishing exponentially and algebraically in the limit the system size $`L\mathrm{}`$. There, the corresponding susceptibility is

$$\chi =L^2m^2/T.$$
(9)
To detect clock symmetry breaking, we consider

$$c_6=\frac{m^6\mathrm{cos}(6\theta )}{m^6}.$$
(10)
It is easy to check that $`c_6`$ equals zero for disordered and KT phases, $`c_6`$ equals $`1`$ for the $`(+0)`$ phase, and $`c_6`$ equals $`1`$ for the $`(+)`$ phase. We have chosen not directly to average $`\mathrm{cos}(6\theta )`$ as its value fluctuates most strongly when $`m`$ is small, that is to say, the ordering we are trying to determine the details of is weakest. The exponent $`m^6`$ has been chosen, somewhat arbitrarily, as the corresponding term in the LGW action is $`\psi _+^6+\psi _{}^6=m^6\mathrm{cos}(6\theta )`$.

### IV.2 Binder cumulant

The appropriate Binder cumulant is:binder

$$U=1\frac{m^4}{3m^2^2}.$$
(11)
The Binder cumulant has a scaling dimension of zero. It thus has the advantage of not requiring fitting an unknown leading exponent. In the limit $`L\mathrm{}`$, the Binder cumulant has the following behaviour: $`U_L0`$ at disordered phase, $`U_L2/3`$ at ordered phase, and $`U_LU^{}`$ at a critical point. For an extended critical phase, the value of $`U^{}`$ can drift, so that its value depends on the precise location within the critical phase.

Since the Binder cumulant has a scaling dimension of zero, the curves for different system sizes at a critical coupling should fall on a line of points where $`U_L=U_L^{}`$ (at least in the region where corrections to scaling are small). As we describe below, this criterion in fact provides an overestimate of the size of the critical phase.

### IV.3 Locating the KT transitions

Here, we briefly describe the finite-size scaling analysis appropriate for KT transitions. This analysis follows that used by Challa and Landau.challalandau  In order to determine the presence of a KT phase, we check if we can fit our simulations to the scaling forms predicted by KT theory. Kosterlitz74  This has the advantage of restricting the number of fitting parameters, which in the most general case (allowing for critical exponents different from the KT ones) would be too large to be practical.

In particular, the correlation length $`\xi `$ and the susceptibility $`\chi `$ behave as

$`\xi \mathrm{exp}(at^{1/2}),`$
$`\chi \xi ^{2\eta },`$
$`m\xi ^{\eta /2}`$ (12)
where $`a`$ is a nonuniversal constant, $`t=(TT_c)/T_c`$ is a reduced temperature, and $`T_c`$ is a critical temperature.

The finite-size scaling form of order parameter and susceptibility is given by

$`m_L=L^bm_0(\xi /L),`$
$`\chi _L=L^c\chi _0(\xi /L),`$ (13)
where $`m_0`$ and $`\chi _0`$ are unknown universal functions and $`b`$ and $`c`$ are constants. It follows from (12) that in the infinite system size limit, we have $`b=\eta /2`$ and $`c=2\eta `$. Therefore, at a critical point, one has:

$`m_LL^{\eta /2},`$ (14)
$`\chi _LL^{2\eta }.`$ (15)
If we have an extended critical phase, these relations should hold over finite temperature range from the upper critical temperature $`T_2`$ to the lower critical temperature $`T_1`$. The critical exponent $`\eta `$ should vary continuously from $`T_2`$ to $`T_1`$. Plotting $`\mathrm{ln}(m_L)`$ or $`\mathrm{ln}(\chi _L)`$ versus $`\mathrm{ln}(L)`$, we can find the critical exponent $`\eta `$ at any point of the extended critical phase.

We can rewrite Eq. (13) as

$`m_LL^b=m_0(L^1\mathrm{exp}(at^{1/2})),`$ (16)
$`\chi _LL^c=\chi _0(L^1\mathrm{exp}(at^{1/2})).`$ (17)
Eq. (16) is valid for $`T<T_1`$ and (17) is valid for $`T>T_2`$. For an appropriately chosen set of parameters $`a,c,T_1`$ the plot of $`\chi _LL^c`$ versus $`L^1\mathrm{exp}(at^{1/2})`$ should collapse onto a universal curve for different system sizes $`L`$. The same should hold for a plot of $`m_LL^b`$ versus $`L^1\mathrm{exp}(at^{1/2})`$. From such a fit, the upper and lower critical temperatures can be determined.

## V Results

In Fig. 1 we present the phase diagram that we find from Monte Carlo simulations. There are three phases: a disordered phase at high temperatures ($`TJ`$) or large magnetic field strengths $`\mathrm{\Gamma }J`$, an extended KT phase at intermediate temperatures, and an $`(+0)`$ ordered phase at low temperatures.

### V.1 The KT phase at $`\mathrm{\Gamma }/J=0.4`$

Our most complete data set was taken for $`\mathrm{\Gamma }/J=0.4`$, for which we describe our analysis in detail in the next few paragraphs. We have chosen this value of $`\mathrm{\Gamma }`$ as a compromise between the following requirements. Firstly, we want to stay well clear of the zero temperature critical points at $`\mathrm{\Gamma }=0`$ and larger $`\mathrm{\Gamma }=\mathrm{\Gamma }_c`$ and concomitant possible crossover phenomena. Secondly, a high transition temperature is needed to allow us to simulate systems with only moderate extent in the imaginary time direction. Thirdly, as our algorithm slows down as the density of domain walls in the imaginary time direction increases, we would like to choose weak quantum fluctuations, that is to say, small values of $`\mathrm{\Gamma }`$. $`\mathrm{\Gamma }/J=0.4`$, located left of centre of the KT dome (Fig. 1) thus appears to be a sensible choice.

First, Fig. 2 shows the behaviour of the order parameter $`m`$ as a function of inverse temperature $`\beta `$ for different system sizes. One can notice easily that the order parameter has strong dependence on the system size in wide range of temperatures – it decreases with increasing system size.

Fig. 3 is a doubly logarithmic plot of the behaviour of the order parameter as a function of the system size. At high temperatures, it extrapolates to zero, whereas it levels out to a constant at low temperatures. For a wide range of temperatures in between, the data points display linear behaviour (without any evidence of logarithmic corrections), beginning at $`\beta J=2.8`$ down to $`\beta J=9.0`$. The exponent $`\eta `$ varies continuously from $`0.323`$ at $`\beta J=2.8`$ to $`0.092`$ at $`\beta J=9.0`$.

This is indeed close to the range expected from KT finite-size scaling (15), although there is an overestimate of the size of the KT phase when compared to the expected range of critical exponents, which is between 1/4 and 1/9. As in the case of the Binder cumulant below, this is a simple consequence of the fact that the difference between the ordered and the critical phase is indiscernible when the size of the system is much less than the correlation length. As the latter grows near the transition, increasingly large finite-size systems in the ordered phase appear to be critical.

To check whether this critical region is indeed delineated by KT transitions, in Fig. 4 we show the data collapse (Eq. 17) for the upper transition. The data scales very well with the following set of parameters $`\beta _2J=3.5,a=5.36,c=1.736`$. Thus we can conclude that the system has a transition between the disordered phase and KT phase at $`\beta _2J=3.5\pm 0.2`$. The critical exponent $`\eta =0.263\pm 0.015`$ at the transition point. This value of the critical exponent $`\eta `$ is close to the theoretical prediction $`1/4`$, and the overestimate appears to be part of a systematic trend discussed below.

In Fig. 5, we show the lower temperature scaled data (see Eq. 16). The data scales quite well (but not perfectly, and over a narrower range than for the upper transition) with the following set of parameters: $`\beta _1J=8.0,a=1.2,c=0.105`$. The error in determining the critical temperature is larger than in the high temperature case. We can conclude that the system has a transition between the KT phase and the ordered phase at $`\beta _1J=8.0\pm 1.0`$. The critical exponent $`\eta =0.105\pm 0.02`$ at the transition point. This value of the critical exponent $`\eta `$ is again close to the theoretical prediction $`1/9`$.

To determine the nature of the ordered phase, we consider the behaviour of $`\mathrm{cos}(6\theta )`$ (Eq. 10). The plot of $`\mathrm{cos}(6\theta )`$ as a function of the inverse temperature $`\beta `$ is shown in Fig. 6. $`\mathrm{cos}(6\theta )`$ goes to $`1`$ as the temperature approaches zero. This implies the existence of the $`(+0)`$-phase at low temperatures for $`\mathrm{\Gamma }/J=0.4`$.

We can also determine the lower transition temperature $`T_1`$ from the criterion of $`\eta =1/9`$ at the lower transition point, which yields $`\beta _1J=7.5\pm 0.5`$. At the transition temperature thus determined, one finds a crossing of $`\mathrm{cos}(6\theta )`$ as a function of $`\beta `$ for different system sizes.oshikawa

Next, we consider the flow diagram of the Binder cumulant, depicted in Fig. 7. In the high (low) temperature phase, its value approaches 0 (2/3) as $`L\mathrm{}`$, whereas it levels off to a value $`U^{}`$ in the KT phase which depends on the strength of the coupling. Challa and Landau challalandau  proposed to use this levelling off as a diagnostic for the KT phase. We do indeed observe such a levelling over a wide range of temperature at $`\mathrm{\Gamma }/J=0.4J`$.

This diagnostic again systematically overestimates the extent of the critical phase. This is evidenced by the uppermost curves in Fig. 7. There is an inflection point at large system size where the value of $`U_L`$, having apparently levelled off, starts increasing again. As one approaches the transition from the ordered phase, this point of inflection wanders to increasingly larger system sizes and hence beyond the scope of the simulations.

### V.2 Other values of $`\mathrm{\Gamma }`$

We now describe what we have found for other values of $`\mathrm{\Gamma }`$. At small $`\mathrm{\Gamma }/J1`$, $`T/J1`$, $`J`$ drops out as an energy scale as the Hilbert space is restricted to the classical ground states, with excited states frozen out by an energy gap of $`O(J)`$. The only dimensionless parameter is then $`\beta \mathrm{\Gamma }`$, so that the lines emanating from the origin of the phase diagram are straight. We estimate their slopes to be $`T_{KT}^l/\mathrm{\Gamma }=0.41\pm 0.05`$ for the lower transition, and for the upper: $`T_{KT}^u/\mathrm{\Gamma }=0.86\pm 0.09`$

The higher temperature data collapse does not always give the critical exponents $`\eta `$ at critical temperatures close to $`1/4`$, especially for large $`\mathrm{\Gamma }`$. For example, we get the following upper set of critical temperatures and critical exponents $`\eta `$: $`\beta _cJ=5.8(6),\eta =0.287(25)`$ at $`\mathrm{\Gamma }/J=0.2`$, $`\beta _cJ=3.5(2),\eta =0.263(15)`$ at $`\mathrm{\Gamma }/J=0.4`$, $`\beta _cJ=2.6(2),\eta =0.230(20)`$ at $`\mathrm{\Gamma }/J=0.8`$, and $`\beta _cJ=3.0(3),\eta =0.203(25)`$ at $`\mathrm{\Gamma }/J=1.2`$. For $`\mathrm{\Gamma }/J=0.2`$, the high-temperature data scales very well over a wide range of temperatures. Thus it is difficult to determine the precise value of the upper transition temperature at this $`\mathrm{\Gamma }`$.

The systematic trend of decreasing $`\eta `$ with increasing $`\mathrm{\Gamma }`$ is in accordance with the fact that $`\eta `$ at the zero temperature transition at large $`\mathrm{\Gamma }`$ is different. For the $`d=3`$ XY universality class, $`\eta `$ is in fact close to 0,zinnjustin  so that the increasing proximity of this fixpoint should be expected to show up in a correction of this kind. One can account for the systematic trend of increasing $`\eta `$ with decreasing $`\mathrm{\Gamma }`$ in the same way. Indeed, the frustrated triangular Ising model has a critical point at $`T=0`$, where $`\eta =1/2`$.stephenson

Probably related to this crossover is the fact that, as $`\mathrm{\Gamma }`$ increases, the plots of the Binder cumulant as a function of system size fail to display the clear flattening visible in Fig. 7 for $`\mathrm{\Gamma }/J=0.4`$ until a system size which is substantially larger; this is displayed for $`\mathrm{\Gamma }/J=0.8`$ in Fig. 8.

The lower transition at larger $`\mathrm{\Gamma }`$ poses a problem for the alogrithm we use. Since we used linked lists to store the location of the domain walls in the temporal direction, the algorithm slows down considerably as $`\mathrm{\Gamma }`$ increases. Together with the increasingly slowly decaying correlations in the space direction as the temperature is lowered, this leads to a considerable increase in the correlation times of our algorithm. As a result, the data we were able to obtain were too noisy to permit a useful data collapse. We can, nonetheless, try to identify the location of the transition by tracking where $`\eta `$ crosses 1/9 (Fig. 1)

We note that the classical six-state clock model, which has a phase diagram very similar to the one we find here, has been studied by many authors.challalandau ; clockmodel  As is the case for us, almost all the authors report that it is much more difficult to get data collapse at the lower KT transition than at the higher. Only in a recent paper, Tomita and Okabe,tomitaokabe  using a ‘probability-changing cluster algorithm’, have found the lower transition easily. Given the frustrated and higher dimensional nature of our spin model, such a cluster algorithm is not available to us.

The nature of the ordered phase at low $`\mathrm{\Gamma }`$ is of the $`(+0)`$ type. There is in principle the possibility of a change of sign of $`v_6`$ for entropic reasons as the couplings are varied.blankschtein ; miyashita  As $`\mathrm{\Gamma }`$ increases, it becomes increasingly hard to determine the sign of $`\mathrm{cos}(6\theta )`$ for the system sizes available to us. The largest $`\mathrm{\Gamma }`$ for which we can confidently state that the lower KT transition is into the $`(+0)`$ phase is $`\mathrm{\Gamma }/J=1.2`$. Although this difficulty may be in part due to a decrease in the strength of $`v_6`$, there is no supporting numerical evidence that the sign of $`v_6`$ ever changes.

### V.3 The quantum critical point

We now turn to the properties of the quantum critical point, i.e. to the zero temperature transition from the clock symmetry broken ordered phase to the paramagnet. We can determine its approximate location to be at $`\mathrm{\Gamma }_c/J=1.65\pm 0.05`$. For large $`\mathrm{\Gamma }`$, the magnet is in a quantum paramagnetic phase with a size-independent susceptibility.

We have not directly attempted to determine the properties of this transition as extensive previous simulations on the classical stacked magnet exist. At the time, there were some suggestions that the observed critical exponents were in fact not those of a $`d=3`$ $`XY`$ modelbunker  and perhaps altered due to the presence of another instability.heino89 ; tricrit

The structure of our phase diagram lends support to the $`d=3`$ $`XY`$ universality class scenario via an independent route. By inducing a dimensional crossover through switching on a finite temperature, we find that the highly nontrivial phase diagram is that predicted by the same Landau theory which gives the $`d=3`$ $`XY`$ universality class. This diagnostic is perhaps more robust than a direct determination of the critical exponents, which can be influenced by corrections to scaling or the proximity of other instabilities.

The shape of the phase boundary near the quantum critical point follows from the knowledge of the critical exponents of the quantum phase transition.subirqptrev  The boundaries of the KT phase near $`\mathrm{\Gamma }_c`$ follow the trajectory $`T_{KT}|\mathrm{\Gamma }\mathrm{\Gamma }_c|^{\nu z}`$, where, for the present case, the dynamical critical exponent $`z=1`$ and $`\nu `$ is close to 2/3.zinnjustin

## VI Conclusion

We have demonstrated that the common action of thermal and quantum fluctuations in the triangular lattice transverse field Ising model generate an interesting fluctuation-driven phase diagram including an extended critical phase bordered by a pair of KT transitions. We have employed several diagnostics for the presence of the KT phase and our results are consistent in considerable detail with what one would expect from an analysis based on a $`d=2+1`$ dimensional $`XY`$ symmetric action with a six-fold clock anisotropy for a temporal dimension of tunable extent. This provides indirect support for the $`XY`$ nature of the quantum critical point. The dynamical properties of the critical phase in close proximity to the quantum critical point are a subject worthy of a study in their own right.

The agreement between numerics and theory, however, is not complete. There are some deviations from the predicted critical exponents, which are probably due to a crossover to the zero temperature critical points. As we have used several diagnostics to determine the nature and location of phases and transitions, we are in a position to compare their relative reliability and find that, for the phase transitions in particular, considering the system size independence of the Binder cumulant and the power-law dependence of the magnetisation on system size systematically overestimate the extent of the critical phase.

In summary, the results we have obtained demonstrate once again that frustrated magnets provide a good starting point for finding unconventional phases and phase diagrams. In this particular case, by using a tunable combination of thermal and quantum fluctuations, we have managed to realise a standard model from statistical mechanics, the $`XY`$ model with sixfold clock anisotropy,Jose77  in dimensions $`d=2`$ and $`d=3`$, in terms of another one, namely the Ising model on the triangular lattice. Given other magnets in this class realise unusual order parameters,ftfim  this approach should provide more opportunities for studying exotic phase diagrams based on simple model spin systems.

## Acknowledgements

We would like to thank D. Huse, A. Karlhede and especially S. Sondhi for many useful suggestions and discussions throughout the course of this work. RM is also grateful to P. Chandra and S. Sondhi for collaboration on closely related work. In addition, we thank A. Bunker, S. Fujiki, T. H. Hansson, J. Lidmar, and M. Oshikawa for useful discussions. This work was in part supported by the Ministère de la Recherche et des Nouvelles Technologies with an ACI grant.


text:
# Chemical evolution of the M82 B fossil starburst

## 1 Introduction

Stellar clusters constitute unique laboratories for studying past and present star formation episodes in galaxies. Their ages and metallicities provide valuable insights into the star formation history as well as the chemical evolution of the host galaxy. Star clusters are especially important for understanding the history of galaxies beyond the Local Group. These galaxies are too far away for their stars to be easily resolved individually, so that most of the information is derived from the integrated properties of star clusters, especially of the brightest, most massive ones.

This field has seen impressive developments during the past decade thanks to the advent of the Hubble Space Telescope (HST) and large ground-based telescopes. One of the main contributions to date of the refurbished HST is the discovery of young clusters with typical globular cluster (GC) masses in galaxies beyond the Local Group. While the formation of GCs was previously thought to be restricted to the early history of the Universe, the discovery of such candidate young GCs shows that the formation of massive stellar clusters, with masses $`m>10^410^5M_{}`$, continues today.

Young and/or intermediate age massive clusters, sometimes referred to as “super star clusters” (SSCs) by virtue of their high luminosities, have been detected mostly in interacting and merging galaxies. The accumulation of such observations (e.g. Holtzman et al. 1996, Zhang, Fall & Whitmore 2001, de Grijs et al. 2001) have highlighted a clear causal link between strong starbursts on the one hand and tidal interactions (e.g. the starburst galaxy M82 which is interacting with its spiral neighbour M81) and/or mergers of galaxies (e.g. the merging pair the “Antennae” NGC 4038/4039) on the other hand.

With the possible exception of the Large Magellanic Cloud, M82 is the nearest galaxy hosting multiple populations of such massive stellar clusters, thus making it a strikingly interesting object. Based on HST (O’Connell et al. 1995) and ground-based imaging (Satyapal et al. 1997), a swarm of young stellar clusters have been unveiled in the central regions of M82. The galactic centre is not the only region where prolific star formation has recently taken place, however. A region extending $``$ 400 - 1000 pc along the major axis northeast of the galactic centre, labelled M82 B (nomenclature from O’Connell & Mangano 1978), has long been suspected to have also gone through a major star formation episode several hundred million years ago. De Grijs et al. (2001) obtained HST multi-passband optical/near infrared observations of this fossil starburst site in which they discovered some 100 star clusters.

In this paper, we perform a detailed study of this cluster system, highlighting its cluster formation history and its chemical evolution. The paper is organized as follows. Section 2 briefly describes the irregular galaxy M82 and the fossil starburst site M82 B. Section 3 explains how we use spectral synthesis models to determine simultaneously the age and the metallicity of the M82 B clusters. Section 4 presents the cluster formation history and the chemical evolution of the M82 B region as well as evidence for the dynamical destruction over time of the clusters. Finally, our conclusions are presented in Sect. 5.

## 2 The irregular galaxy M82

M82 is the prototype starburst galaxy, being the nearest and, therefore, the best studied example of this class of galaxies. Freedman et al. (1994) estimated its distance to be 3.6 Mpc (i.e. its distance modulus is $`mM=27.8`$). Such a proximity makes it uniquely valuable for a wide range of studies regarding extragalactic star clusters. M82 belongs to the nearest group of interacting galaxies, whose other prominent members are the large spiral M81 and a fainter elliptical companion, NGC 3077. From optical images only, this set of three galaxies appears quiescent, showing little sign of tidal interactions. However, a completely different conclusion is reached from radio surveys. A map of neutral hydrogen obtained with the Very Large Array (Yun et al. 1994) reveals huge streams of neutral hydrogen gas connecting the three galaxies as well as several fainter ones. These HI streams are dominated by filamentary structures, loops and twists thus suggesting that the three largest members of the M81 group have gone through several close encounters. The corresponding tidal interactions have most probably triggered the formation of some of the M82 star clusters. For instance, they may be responsible for channeling large amounts of gas into the central regions of the galaxy, thus inducing the active starburst in the central regions (Rieke et al. 1993). The study of the off-centre region M82 B is particularly valuable since it has been shown to contain a more evolved star cluster system, and to have the characteristics of an older burst of star formation (e.g., O’Connell & Mangano 1978, Marcum & O’Connell 1996, de Grijs et al. 2001). Therefore, the study of M82 B will enable us to probe the galaxy’s history to earlier times.

Starburst regions are often obscured so heavily by dust that their study at infrared wavelengths represents a definitive advantage with respect to the optical. De Grijs et al. (2001) used the Wide Field Planetary Camera 2 (WFPC2) as well as the Near-Infrared Camera and Multi-Object Spectrometer (NICMOS), on board HST, to obtain high resolution imaging of M82 B in both the optical and the near infrared. Two adjacent 35” square fields, named “B1” and “B2”, were imaged, the “B2” region being the closest to the galactic centre and the active starburst. Despite the fact that B1 and B2 simply correspond to two adjacent fields, this differentiation, i.e. M82 B1 and M82 B2 instead of M82 B, will often be used throughout this paper as both regions show distinct differences in their chemical evolution (see Sect. 4.4). Their HST observing run provided de Grijs et al. (2001) with the $`B`$, $`V`$, $`I`$ and $`J`$ magnitudes of the star clusters detected in M82 B. The availability of these four magnitudes allows us to disentangle the extinction, age and metallicity effects. While the way of estimating simultaneously the age and the metallicity is presented in Sect. 3.1, we now briefly summarize how de Grijs et al.(2001) derived the intrinsic colour indices of the M82 B clusters, these colour indices being used in our study.

In a $`(BV)`$ vs $`(VI)`$ diagram, the aging trajectories (i.e. the colour evolution with time at a given metallicity) and the extinction vector are not entirely degenerate for this choice of optical colors, showing that the $`BVI`$ diagram could be used to estimate the colour excesses by which a cluster is reddened (de Grijs et al. 2001, their Fig. 12). This method also allows to derive extinction values for each cluster individually. This is useful in a galaxy like M82, which has a significant and highly variable extinction and which cannot be characterized by a single mean value for the internal extinction and reddening. Assuming that the metallicity of the M82 B clusters is roughly solar, de Grijs et al. (2001) used this property to estimate their intrinsic colour indices, their extinction and their absolute magnitudes.

Using the colour indices dereddened as explained above, de Grijs et al. (2001) determined the age distribution of the M82 B clusters, which contains information about the cluster formation rate. The age of the M82 B clusters ranges from 200 Myr to over 10 Gyr, a time interval over which clusters have been forming continuously, although not at a constant rate. The most striking feature is a strong peak of cluster formation starting about 1 Gyr ago and lasting about 600 Myr (de Grijs et al. 2001), a time roughly corresponding to the last perigalactic passage of M81 and M82 (Brouillet et al. 1991). The close match between both of these independently determined time-scales provides, therefore, very strong evidence that interactions between galaxies lead to enhanced star cluster formation.

While de Grijs, Bastian & Lamers (2002) study the dynamical history of the M82 B star cluster system in great detail, we focus on its chemical evolution. In order to carry out such an analysis, we refine the results obtained by de Grijs et al. (2001), that is, we redetermine both the age and the metallicity of the star clusters.

## 3 Spectral Synthesis

Since the pioneering papers of Tinsley (1972) and Searle, Sargent & Bagnuolo (1973), spectral population synthesis, i.e. the modeling of the spectral energy distribution of a star cluster, has become a standard technique to study the integrated stellar populations of galaxies. Assuming that stellar clusters are single-age and single abundance groups of stars, their integrated colours reflect their age and metallicity (for a given initial mass function, IMF).

In order to perform a self-consistent comparison of our results with those obtained by de Grijs et al. (2001), we use the same spectral synthesis model, i.e. the Bruzual & Charlot (1996, hereafter BC96) model. We consider the photometric properties of Simple Stellar Populations (SSPs; i.e. single-age and single abundance groups of stars) with a Salpeter IMF from 0.1 M up to 125 M. BC96 provides the temporal evolution of the broad-band colours of such SSPs for a wide range of metallicities. In the following, we consider the three intrinsic colour indices $`(BV)_0`$, $`(VI)_0`$ and $`(VJ)_0`$ whose observational counterparts were obtained by de Grijs et al. (2001). We show below how the addition of the $`(VJ)_0`$ colour index enables us to estimate both the age and the metallicity of the clusters.

### 3.1 Disentangling age and metallicity effects

A greater age or a higher metallicity both lead to redder colours, an effect known as the age-metallicity degeneracy. For instance, in a $`(BV)_0`$ versus $`(VI)_0`$ plot, the iso-metallicity tracks (i.e. the locus of an SSP at different ages for a given metallicity) and the isochrones are parallel and located in the same part of the optical colour-colour diagram, precluding therefore any simultaneous estimate of age and metallicity. De Grijs et al. (2001) also obtained the J-band magnitude for each cluster. Figure 1 explores the suitability of a $`VIJ`$ colour-colour diagram to lift this age-metallicity degeneracy. Thin lines represent isometallicity tracks, i.e. aging trajectories, and thick lines refer to isochrones (from left to right, log $`t/\mathrm{yr}`$ = 8.5, 9.0, 9.5 and 10.0). The grid defined by the isometallicity tracks and the isochrones shows that a $`(VJ)_0`$ vs $`(VI)_0`$ diagram is well suited to disentangle the effects of age and metallicity.

As mentioned in Sect. 2, the cluster intrinsic colours were determined by de Grijs et al. (2001) based on the hypothesis that the metallicity of the clusters in M82 B was roughly solar. However, the next section shows that M82 B hosts very metal-poor (Z $``$ 0.02Z) and super-solar abundance (Z $`>`$ Z) clusters in addition to solar abundance ones. Therefore, one can question the validity of the present analysis which makes use of the intrinsic colours determined by de Grijs et al. (2001). The solar metallicity track and the extinction trajectories being not entirely degenerate in a $`(BV)_0`$ vs $`(VI)_0`$ diagram, de Grijs et al. (2001) were able to estimate the visual extinction and the colour excess for each cluster (see Sect. 2). Furthermore, the iso-metallicity tracks are at the same locus of the diagram, with the distance between two subsequent tracks of distinct metallicity being smaller than the photometric uncertainties. As a consequence, the visual extinction values and the colour excesses (or, equivalently, the intrinsic colours) estimated by de Grijs et al. (2001) on the basis of the solar metallicity hypothesis, can be used throughout this work irrespective of the metallicity of the clusters.

### 3.2 Accuracy of the age and metallicity estimates for unextinguished clusters

The medians of the $`1\sigma `$ error distributions for the different colour indices are $`\sigma _{(BV)}`$ 0.08, $`\sigma _{(VI)}`$ 0.08 and $`\sigma _{(VJ)}`$ 0.15. Let us first assume that the uncertainties of the observed colours constitute a good approximation of the intrinsic colour uncertainties. The typical uncertainties in $`(VI)`$ and $`(VJ)`$ are superimposed on the age-metallicity grid in Fig. 1 for the solar metallicity SSPs of 600 Myr and 10 Gyr. The corresponding uncertainties in age are the smallest for SSPs with ages from 300 Myr to 1 Gyr. For instance, for an SSP of solar metallicity the typical error in the estimated age is not greater than 40% at an age of 600 Myr, while it is as large as a factor of 2 for a 10 Gyr old SSP. This greater uncertainty at greater age, which is related to the tightening of the isochrones at log $`t/\mathrm{yr}9.5`$ in the $`VIJ`$ diagram, is another aspect of the decrease in sensitivity of the colour indices to ages as an SSP ages (e.g. Worthey et al. 1994). Interestingly, the diagram shows the greatest sensitivity to ages within the time interval of interest to us, namely the age of the burst of cluster formation at $`\mathrm{log}t/\mathrm{yr}9`$.
As for the accuracy of the metallicity estimates, the differences in the $`(VJ)_0`$ colour index between the isometallicity tracks Z = Z, Z = 0.4Z and Z = 0.2Z are of a similar order of magnitude as the uncertainty $`\sigma _{(VJ)}`$ so that we cannot easily distinguish among such close-to-solar metallicities based on the available data set. However, the $`VIJ`$ colour-colour diagram can, in fact, be used to distinguish solar metallicity SSCs from very metal-poor (i.e. Z = 0.02Z) and super-solar abundance (i.e. Z $`>`$ Z) populations.

We caution however that this reasoning does not take into account the extinction uncertainties. Figure 1 displays the extinction vector in the $`VIJ`$ diagram and clearly shows that extinction uncertainties cannot be neglected as they will combine with the colour uncertainties quoted above to enlarge the error boxes in Fig. 1. This reasoning therefore applies only to clusters whose intrinsic and observed colours are not significantly different. In Sect. 4.2, we will estimate in a more careful way how the magnitude uncertainties (i.e. $`\sigma _B`$, $`\sigma _V`$, $`\sigma _I`$, $`\sigma _J`$) propagate as extinction, age, mass and metallicity uncertainties.

## 4 The Cluster Formation History of M82 B

In this section, we estimate the age, the metallicity and the mass-to-light ratio of the M82 B clusters by comparing their extinction-corrected colours to those predicted as a function of age and metallicity by BC96. We will then derive the cluster formation history of the M82 B region.

### 4.1 Age, metallicity and mass estimates

The age ($`t_{\mathrm{cl}}`$) and the metallicity ($`Z_{\mathrm{cl}}`$) of each stellar cluster are estimated simultaneously by means of a least-squares minimisation:

$$\mathrm{min}\left\{S^2(t,Z)=\underset{i=1}{\overset{3}{}}\left(\frac{\mathrm{CI}_{i}^{}{}_{}{}^{\mathrm{intr}}\mathrm{CI}_{i}^{}{}_{}{}^{\mathrm{SSP}}(t,Z)}{\sigma _{\mathrm{CI}_i}}\right)^2\right\}$$
(1)
$``$ tcl and Zcl .

In this relation, $`\mathrm{CI}_{i=1,2,3}^{\mathrm{intr}}`$ denote the three intrinsic colour indices (i.e. $`(BV)_0`$, $`(VI)_0`$, $`(VJ)_0`$) and CI$`{}_{i=1,2,3}{}^{\mathrm{SSP}}(t,Z)`$ represent the corresponding BC96 integrated colours of an SSP with age $`t`$ and metallicity $`Z`$. $`\sigma _{\mathrm{CI}_{i=1,2,3}}`$ are the 1$`\sigma `$ uncertainties in the colour indices. The resulting age and metallicity distributions are shown in Figs. 2 and 3, respectively. They will be discussed, along with the other results, in the next sections. The different panels of Fig. 4 display the comparison between the BC96 colour indices for the derived ages and metallicities and the corresponding intrinsic colour indices obtained from the multicolour HST photometry. The 1$`\sigma `$ error bars and, when relevant, the upper/lower limits (indicated by the arrows) of the observed colours are included. The good agreement between the observed colour indices and those derived from the BC96 spectral synthesis models illustrates the reliability of our age and metallicity estimates.

Using these age and metallicity estimates, we can now apply the age and metallicity dependent mass-to-light ratio $`\frac{M}{L_V}(t,Z)`$ predicted by BC96 for a single burst stellar population to derive masses $`m`$ for our cluster sample.

As quoted above (Sect. 3), the mass-to-light ratio predicted by the BC96 model we use assumes a Salpeter IMF. However, there is now ample evidence that a single power-law Salpeter IMF can be ruled out at masses lower than 1 M, the IMF slope becoming clearly much shallower than the Salpeter IMF. For instance, the IMF observed in the solar neighbourhood may be roughly flat at the lower end of the stellar mass spectrum or may even peak at a mass around 0.25 M and decline into the brown dwarf regime (Larson 1998). The Salpeter IMF obviously overestimates the relative number of subsolar-mass stars. These stars being those with the largest mass-to-light ratio, a Salpeter IMF causes an overestimate of the individual cluster masses, although the relative mass distribution of our entire cluster sample remains unaffected. If we had used more modern IMFs (e.g. Kroupa, Tout & Gilmore 1993) this would have resulted in lower cluster mass estimates, by a factor $``$ 2-3 (see de Grijs et al. 2002) with respect to the mass estimates presented in this paper.

Our best age, mass and metallicity estimates for each cluster, along with their uncertainties (see the following section) are listed in Tables 1 and 2 for B1 and B2, respectively. The results are presented in Figs. 2 to 9. In these figures, the few clusters for which either $`B`$, $`V`$ or $`I`$ is an upper or lower limit only are represented by open symbols or dashed histograms.

### 4.2 Uncertainty estimates

We now address the issue of estimating the errors by which the derived cluster properties are affected. Since the extinction was derived from the $`BVI`$ diagram, it is clear that uncertainties in the optical filters will propagate as extinction uncertainties. In the $`BVI`$ colour-colour diagram, the extinction is the distance, measured along the reddening trajectory, from the observed point ($`(BV)`$,$`(VI)`$) down to the aging trajectories (Fig. 12 of de Grijs et al. 2001). Considering stellar populations older than 100 Myr, the age range of the vast majority of the clusters of our sample, the aging trajectories display the same linear growth of $`(VI)_0`$ vs $`(BV)_0`$, whatever the metallicity. It is therefore reasonable to approximate all the iso-metallicity tracks by the same straight line with slope and ordinate origin $`m_{\mathrm{SSP}}`$ ($`0.87`$) and $`p_{\mathrm{SSP}}`$ ($``$ 0.28), respectively:

$$(VI)_0m_{\mathrm{SSP}}(BV)_0+p_{\mathrm{SSP}}.$$
(2)
Under this assumption, it is easy to show that the extinction modulus obeys:

$$\begin{array}{c}A_V\frac{\sqrt{1+m_{\mathrm{ext}}^2}}{\sqrt{E(BV)^2+E(VI)^2}}\hfill \\ \hfill \times \frac{(VI)m_{\mathrm{SSP}}(BV)p_{\mathrm{SSP}}}{m_{\mathrm{ext}}m_{\mathrm{SSP}}}.\end{array}$$
(3)
In this equation, $`E(BV)`$=0.32 and $`E(VI)`$=0.52 are the optical colour excesses corresponding to an extinction of one magnitude, $`m_{\mathrm{ext}}`$ is the slope of the reddening trajectory, i.e. $`m_{\mathrm{ext}}=0.52/0.32`$ (Cardelli, Clayton & Mathis 1989). Substituting these numerical values into Eq. (3), we get an estimate of the extinction dimming a cluster whose measured optical magnitudes are $`B`$, $`V`$ and $`I`$:

$$A_V4.14\left[(1+m_{\mathrm{SSP}})Vm_{\mathrm{SSP}}BIp_{\mathrm{SSP}}\right].$$
(4)
If we assume that the $`BVI`$ measurements are entirely independent, the $`1\sigma `$ error by which the extinction is affected derives from Eq. (4):

$$\sigma _{A_V}4.14\sqrt{m_{\mathrm{SSP}}^2\sigma _B^2+(1+m_{\mathrm{SSP}})^2\sigma _V^2+\sigma _I^2}.$$
(5)
For each cluster, Tables 1 and 2 also list the best, minimal and maximal values of the extinction $`A_v`$. The observed points of a few clusters fall below the aging trajectories in the $`BVI`$ diagram but remain nevertheless consistent with the theoretical tracks within the photometric uncertainties (de Grijs et al. 2001, their Fig. 12). In such cases, we have assumed zero extinction.

Figure 1 displays the reddening vector corresponding to an extinction of one magnitude in the $`VIJ`$ diagram, i.e. the diagram used to disentangle the age and the metallicity. Since this vector is roughly parallel to the iso-metallicity tracks, the extinction uncertainties (Eq. (5)) will significantly propagate as age uncertainties. An overestimate of the extinction leads to an underestimate of the cluster age, and vice versa. In addition, the age dependent mass-to-light ratio causes the age uncertainties to propagate as mass uncertainties.

In order to perform a consistent estimate of the age uncertainties, we propagate the photometric uncertainties in each filter $`B`$, $`V`$, $`I`$ and $`J`$. The age $`t`$ of a cluster is a function of its $`BVIJ`$ photometry, that is, $`(\mathrm{log}t)=\mathrm{log}t(B,V,I,J)`$. Assuming that $`B`$,$`V`$,$`I`$ and $`J`$ are independent measurements, $`\sigma _{\mathrm{log}t}`$, the $`1\sigma `$ error in $`\mathrm{log}t`$ we aim to estimate, is given by:

$$\begin{array}{c}\sigma _{\mathrm{log}t}^2=\sigma _B^2\left(\frac{(\mathrm{log}t)}{B}\right)_{V,I,J}^2+\sigma _V^2\left(\frac{(\mathrm{log}t)}{V}\right)_{B,I,J}^2\hfill \\ \hfill +\sigma _I^2\left(\frac{(\mathrm{log}t)}{I}\right)_{B,V,J}^2+\sigma _J^2\left(\frac{(\mathrm{log}t)}{J}\right)_{B,V,I}^2.\end{array}$$
(6)
In this equation, $`\sigma _B`$ is the $`B`$ magnitude uncertainty, etc, $`((\mathrm{log}t)/B)_{V,I,J}`$ is the partial derivative of the age logarithm with respect to $`B`$ while keeping $`V`$,$`I`$,$`J`$ constant, etc. Equation (6) can be rewritten as:

$$\begin{array}{c}\sigma _{\mathrm{log}t}^2=[\delta _B(\mathrm{log}t)]^2+[\delta _V(\mathrm{log}t)]^2\hfill \\ \hfill +[\delta _I(\mathrm{log}t)]^2+[\delta _J(\mathrm{log}t)]^2\end{array}$$
(7)
where $`\delta _B(\mathrm{log}t)`$ is the variation in $`\mathrm{log}t`$ caused by a 1$`\sigma `$ change in the $`B`$ magnitude while keeping the $`V`$, $`I`$, $`J`$ ones to their best values, etc. The different terms of the right-hand side of Eq. (7) can therefore be derived by computing the cluster age for each of the following combinations of observed magnitudes: $`(\mathrm{log}t)(B,V,I,J)`$, $`(\mathrm{log}t)(B\pm \sigma _B,V,I,J)`$, $`(\mathrm{log}t)(B,V\pm \sigma _V,I,J)`$, $`(\mathrm{log}t)(B,V,I\pm \sigma _I,J)`$, $`(\mathrm{log}t)(B,V,I,J\pm \sigma _J)`$, $`B`$, $`V`$, $`I`$ and $`J`$ being the best magnitude values. Figure 5 shows how such 1$`\sigma `$ individual changes in the observed magnitudes of a cluster affect its location in the $`VIJ`$ diagram. The left and right panels illustrate the cases of B1-43 and B2-32, respectively. The crossed-circle indicates the location of the cluster corresponding to its best magnitude values. The diagrams illustrate that changing these values by the following amounts:

* $`+\sigma _B`$ or $`\sigma _V`$ or $`+\sigma _I`$ or $`+\sigma _J`$ (corresponding to the stars in the panels of Fig. 5) leads to an increase of the age estimate. Equation (4) indeed indicates that increasing $`B`$ or $`I`$, or decreasing $`V`$ leads to a decrease of the extinction and, therefore, to a shift towards the right part of the $`VIJ`$ diagram, i.e. towards older ages. Increasing the $`J`$ magnitude does not affect the extinction but leads nevertheless to an age increase through a downwards shift in the $`VIJ`$ diagram;
* a similar reasoning shows that changing the magnitudes by $`\sigma _B`$ or $`+\sigma _V`$ or $`\sigma _I`$ or $`\sigma _J`$ (the diamonds in the panels of Fig. 5) with respect of their best values lowers the cluster age estimate.

Summing quadratically the age changes derived from the series (a) will therefore provide the age “upper” error bar, while summing quadratically the age changes derived from the series (b) will provide the “lower” age error bar, that is:

$$\begin{array}{c}(+\sigma _{\mathrm{log}t})^2=[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B+\sigma _B,V,I,J)]^2\hfill \\ \hfill +[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B,V\sigma _V,I,J)]^2\\ \hfill +[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B,V,I+\sigma _I,J)]^2\\ \hfill +[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B,V,I,J+\sigma _J)]^2\end{array}$$
(8)
and

$$\begin{array}{c}(\sigma _{\mathrm{log}t})^2=[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B\sigma _B,V,I,J)]^2\hfill \\ \hfill +[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B,V+\sigma _V,I,J)]^2\\ \hfill +[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B,V,I\sigma _I,J)]^2\\ \hfill +[(\mathrm{log}t)(B,V,I,J)(\mathrm{log}t)(B,V,I,J\sigma _J)]^2\end{array}$$
(9)
for the “upper” and “lower” error bars in $`\mathrm{log}t`$, respectively.

As already mentioned, the age uncertainties propagate as mass uncertainties through the age dependent mass-to-light ratio. The errors affecting the cluster mass and metallicity are derived in a way similar to those affecting the age: we individually propagate the magnitude uncertainties and we then compute the lower/upper error bars by summing quadratically the shifts leading to mass and metallicity estimates smaller/greater than the best ones. The 1$`\sigma `$ mass error we derive thus takes into account the uncertainties affecting both the intrinsic luminosity of the cluster and its mass-to-light ratio. The luminosity uncertainties originate from

1. the visual magnitude uncertainty ($`\sigma _V`$),
2. the extinction uncertainty, which depends on the optical magnitude uncertainties, i.e. $`\sigma _B`$, $`\sigma _V`$ and $`\sigma _I`$ (see Eq. (5)).

The mass-to-light ratio uncertainties are driven by the age and metallicity uncertainties, which both depend on the $`BVIJ`$ photometry. The metallicity uncertainties however mostly reflect the $`J`$ magnitude 1$`\sigma `$ error. More than two thirds of our clusters are characterized by:

$$\begin{array}{c}Z_{\mathrm{max}}=Z(B,V,I,J\sigma _J)\mathrm{and}\hfill \\ \hfill Z_{\mathrm{min}}=Z(B,V,I,J+\sigma _J).\end{array}$$
(10)
This result is not unexpected: the extinction vector being roughly parallel to the aging trajectories, the influence of the extinction uncertainties (driven by the optical uncertainties) upon the metallicity is strongly reduced.

The age, mass and metallicity estimates (minimal, best and maximal values) are listed in Tables 1 and 2 for B1 and B2, respectively. Our photometric mass estimates show an accuracy of well within an order of magnitude. Tables 1 and 2 show that the mass estimates are more robust than the age estimates (see also Fig. 8). This pattern comes from the combination of two effects. Firstly, the mass-to-light ratio grows with time at a slower rate than time itself (i.e. the evolution of $`\mathrm{log}(m/L_v)`$ vs $`\mathrm{log}t`$ exhibit a linear behaviour with a slope of less than unity, see BC96). In addition, the impact upon the mass estimate of mass-to-light ratio variations is partially counterbalanced by extinction variations. In order to illustrate this, let us consider a cluster for which the extinction has been overestimated. The adjustment of the extinction to a more appropriate (i.e. lower) value leads to a smaller cluster intrinsic brightness which, alone, would lower the photometric mass estimate. At the same time however, the extinction decrease causes a rise of the cluster age estimate (see the extinction vector in Fig. 5) and, therefore, an increase of the cluster mass-to-light ratio. These effects, a luminosity decrease and a mass-to-light ratio increase, act upon the mass estimate in opposite ways. For instance, considering a roughly solar metallicity model, a decrease of 0.6 magnitude in $`A_V`$ causes an age increase of $`\mathrm{log}t0.5`$ at the time of the burst (see Fig. 5). Accordingly, the mass-to-light ratio grows by a factor of $``$ 2.7. The new mass estimate is thus larger than the previous one by a factor 2.7/2.5$`{}_{}{}^{0.6}1.6`$. In summary, while an age increase leads to a growth of the mass estimate through the age dependent mass-to-light ratio, this mass increase is reduced by a decrease in extinction, which acts in the opposite way. It is worth to mention that this extinction-induced reduction of the mass variation is not so important for intermediate and old populations. In fact, due to the tightening of the isochrones at older age (see Fig. 5), a given extinction variation corresponds to a greater age range at older age than at the burst epoch and, therefore, to a larger mass-to-light ratio variation. The relative reduction of the mass uncertainty thanks to the extinction variation is thus much weaker at old age than at the burst age.

As a result of the variations with time of the mass-to-light ratio, the lower and upper mass limits are strongly coupled to the lower and upper age limits, respectively.

### 4.3 The cluster formation history

Figure 2 shows the age distributions of the clusters in the entire M82 B region (top panel), in regions B1 (middle panel) and B2 (bottom panel). Clusters have been forming continuously, their age ranging from more than 10 Gyr down to about 10 Myr. In addition, the B1 and B2 regions exhibit a period of enhanced cluster formation, from about 1.5 Gyr ago until 500 Myr ago. We therefore confirm the main result of de Grijs et al. (2001), although our peak of cluster formation is slightly older, $``$ 1 Gyr (see also de Grijs et al. 2002, for more details). The M82 B burst was followed by a sharp decline in the cluster formation rate. Less than one fifth of the clusters are younger than 500 Myr. Starbursts are likely to be strongly self-limited by supernova-driven outflows which remove the remaining cool gas from the immediate starburst region. It is interesting to note that the B2 region, which is closer to the active starburst, has formed a few clusters more recently than B1.

It has long been suspected that the tidal interactions between the M81 group members (i.e. M81, M82 and NGC 3077) are responsible for the starburst nature of M82 and for the HI bridges connecting these galaxies. Brouillet et al. (1991) used an N-body model, representing M81, M82 and NGC 3077, to simulate the deformations of the neutral hydrogen distribution due to the tidal interactions: their results reproduce the observed HI bridges and tails rather well. According to their simulation, the last perigalactic passage of M82 and M81 took place some 500 Myr ago. The onset of the burst about 1.5 Gyr ago is thought, therefore, to coincide with the beginning of the tidal interactions between M82 and its prominent neighbour M81.

### 4.4 The chemical evolution

The M82 B cluster system is characterized by a wide range of metallicities (see Fig. 3) as it hosts both very metal-poor clusters (Z$``$0.02Z) and super-solar abundance ones (Z$`>`$Z). Such a pattern must be related to the region’s chemical evolution. In order to illustrate this, each panel of Fig. 6 shows the cluster age distribution at a given metallicity for regions B1 (left column) and B2 (right column). The metallicities shown are, from top to bottom, $`Z=0.02Z_{}`$, $`0.2Z_{}`$, $`0.4Z_{}`$, $`1Z_{}`$, $`2.5Z_{}`$, $`5Z_{}`$.

From the birth of its first clusters until about 1 Gyr ago, M82 B formed clusters with roughly solar and subsolar metal abundances only. In both subregions B1 and B2, the formation of the first clusters with significantly super-solar (i.e. $`Z>Z_{}`$) metallicities coincides with the onset of the burst. The stellar ejecta from the evolution of older stars have most probably contributed to the chemical enrichment.
While this chemical evolution up to super solar metallicities at the burst epoch is a feature in common for B1 and B2, both subregions nevertheless exhibit a striking difference regarding the most metal-poor clusters. At the time of the burst, and despite the chemical enrichment noticed above, B2 managed to form a significant subpopulation of very metal-poor clusters (i.e. $`Z0.02Z_{}`$) with respect to B1 (top panels in Fig. 6). We stress that half of these B2 clusters have very well-constrained ages ($`(\mathrm{log}t_{\mathrm{max}}\mathrm{log}t_{\mathrm{min}})0.5`$) and metallicities ($`Z_{\mathrm{max}}=0.2`$), so that we probably do not face a spurious effect.

Either some interstellar gas in B2 managed to escape the chemical enrichment driven by the older stellar populations, or alternatively some “fresh” metal-poor gas was injected into the B2 region shortly before the burst. Whereas the first scenario is hard to explain, the second is more appealing. A large amount of circumgalactic cold material is orbiting M82, stripped from the gas-rich outskirts of the M81 group spirals by tidal interactions. These interactions involve quite a few galaxies, of which three are of a non-negligible size (i.e. M81, M82 and NGC 3077). It is therefore not unreasonable to assume that tidal interactions and gas tidal stripping were at least somewhat important even before M81 and M82 almost collided about 0.5 Gyr ago. If the bridges of circumgalactic gas were indeed already present 1 Gyr ago (i.e. 500 Myr prior to the perigalactic passage), gas infall from this reservoir onto the B2 region may have induced the formation of this subpopulation of metal-poor clusters. Alternatively, some metal-poor gas associated with the M82 outskirts could have been driven onto B2 by tidal interactions.

In comparing the properties of the clusters of B1 with those of B2, we implicitely assume that their respective current and birth locations, with respect to the galactic centre, coincide. The validity of this assumption may be questioned as, at the distance from the centre of region B, one would expect M82’s differential rotation (Shen & Lo 1995) to have caused the starburst area to disperse. However, the lack of bright clusters like those in M82 B outside this region does not meet this expectation and, on the contrary, indicates that the fossil starburst region has remained relatively well constrained. The reason for this is likely found in the complex structure of the disc. It is well-known that the inner $`1`$ kpc of M82 (i.e. the radial extent of M82 B) is dominated by a stellar bar (e.g., Wills et al. 2000) in solid-body rotation. From observations in other galaxies, it appears to be a common feature that central bars are often surrounded by a ring-like structure. If this is also true for M82, it is reasonable to assume that stars in the ring are trapped, and therefore cannot move very much in radius because of dynamical resonance effects. As a consequence, not only has the fossil starburst site managed to stay together, both subregions B1 and B2 may also have conserved their own specificity for about 1 Gyr (see also de Grijs 2001). In the next section, we show that the presence (B2), or the absence (B1), of a subpopulation of very metal-poor clusters at the burst epoch is not the only difference between B1 and B2.

#### 4.4.1 Did the B1 clusters go through a self-enrichment phase?

A high-pressure medium favours the formation of pressure-bound clouds, which may be the progenitors of future stellar clusters. For instance, Fall & Rees (1985) suggested that the formation of the halo GCs in our Galaxy occured in dense and cold clouds which are in pressure equilibrium with a hot and diffuse background. This two-phase medium is thought to have been formed by the collapse of the protoGalaxy. In the frame of their theory, the GC progenitor clouds are thermally supported and made of primordial gas. To explain the metallicities of the Galactic halo GCs, Parmentier et al. (1999) further developed this picture for GC formation and proposed that the halo GC gaseous progenitors went through a self-enrichment phase. A self-enrichment scenario assumes the formation of a first generation of stars in the central regions of each proto-globular cluster cloud. When the massive stars explode as Type II supernovae, they chemically enrich the surrounding cloud gas and trigger the formation of a second stellar generation, more metal-rich than the first one. These second generation stars form the proto-globular cluster. Supernovae having long been thought to disrupt the cloud of gas out of which they have formed, Parmentier et al. (1999) studied the ability of pressure-truncated clouds to retain supernova ejecta and further showed that such a process is able to explain Galactic halo GC metallicities. Assuming that the number of supernovae is the maximum number predicted by their model (i.e. a few hundred), a self-enrichment episode in pressure-bound clouds shows up as a correlation between the mass and the achieved metallicity of the progenitor cloud in the sense that the least massive clouds are the most metal-rich. Such a trend emerges because if the bound pressure is higher, the mass of the pressure-truncated cloud will be lower, and its ability to retain supernova ejecta will be greater. Regarding the resulting system of GCs, a tight mass-metallicity correlation should not be expected, however. The star formation efficiency (i.e. the ratio between the mass of the second stellar generation and the mass of gas) and the number of supernovae (which determines the amount of metals dispersed within the cloud) may vary from cloud to cloud, imprinting therefore a scatter on the initial mass-metallicity correlation (see Parmentier & Gilmore 2001 for a more detailed discussion). In spite of this scatter, Parmentier & Gilmore (2001) showed that such a mass-metallicity correlation is statistically present in the Milky Way Old Halo, i.e. the group of old and coeval halo GCs (e.g. Rosenberg et al. 1999). At this stage, it is worth keeping in mind that the search for a mass-metallicity relation must be restricted to a coeval population of clusters. The (surviving) correlation primarily reflects the variations in the pressure by which the clouds are bound and has nothing to do with the chemical enrichment with time.

Both the M82 interstellar medium and the hot protogalactic background are high-pressure mediums. Therefore, a tempting question is whether the almost coeval population of clusters formed during the burst (0.5 $`<t<`$ 1.5 Gyr) have undergone a self-enrichment process in pressure-truncated progenitor clouds. These M82 B clouds were most likely mainly supported by non-thermal effects (e.g. turbulence and magnetic fields) and resembled more closely giant molecular clouds than the neutral hydrogen clouds described by the Fall & Rees (1985) theory. As mentioned above, however, the key parameter leading to a trend between the mass and the metallicity of stellar clusters is the pressure of the medium in which their progenitor clouds were embedded. It is widely believed that the boundary of giant molecular clouds is set by pressure balance with the surrounding, more diffuse, interstellar medium (e.g. McLaughlin & Pudritz 1996). Therefore, the hypothesis of a mass-metallicity relation among the clusters of the M82 B system makes sense. In order to explore the possibility of a self-enrichment event (superimposed on an unavoidable pre-enrichment phase owing to the young age of the clusters), we now search for such a mass-metallicity correlation among the $``$1 Gyr old clusters in the B1 and B2 regions. Again, B1 and B2 exhibit different patterns (see Fig. 7, B1/B2: top/bottom panel). The B1 sample shows a mass-metallicity correlation in the sense expected by the self-enrichment model, i.e. the most metal-rich clusters are, on average, least massive: the linear Pearson correlation coefficient is $`r=`$0.70, corresponding to a correlation probability of $`\mathrm{}`$ = 99.5%. The scatter in metallicity at a given mass interval (e.g. $`10^5<m<3\times 10^5M_{}`$) remains significant but this is not unexpected because of the scatter sources mentioned above. The age of the vast majority of the B1 burst clusters is well constrained ($`(\mathrm{log}t_{\mathrm{max}}\mathrm{log}t_{\mathrm{min}})0.7`$ for 11 of the 13 clusters) and so is their mass (see the horizontal error bars in the top panel of Fig. 7). In addition, the metallicity error bars are rather limited as well. As a result, the error bars, both in mass and metallicity, do not significantly corrupt the correlation.
In contrast to B1, the B2 burst clusters do not exhibit any mass-metallicity correlation. Removing the most metal-poor clusters (i.e. those assumed to have formed from circumgalactic gas) does not change the situation. This difference between the B1 and B2 mass-metallicity diagrams might be related to the formation of some metal-poor clusters in B2 at the burst epoch. If these low metal abundance clusters indeed formed out of the arrival of “fresh” gas, a shock wave would have accompanied the gas injection, perturbing the external pressure and preventing therefore any self-enrichment in pressure-bound clouds.
While the scenario presented above is tentative, the fact that B1 and B2 exhibit different behaviours from the points of view of the chemical evolution (top panels of Fig. 6) and of the mass-metallicity relation (Fig. 7) is puzzling and suggests actual differences between the cluster formation histories in both regions at the time of the burst.

### 4.5 Dynamical disruption of the SSCs

Figure 8 shows the evolution with time of the mass of the clusters in B1 and B2. The most striking feature is the deficit in low-mass clusters (i.e. $`\mathrm{log}m/M_{}<5.5`$) during the pre-burst phase (log $`t`$/yr $`9.2`$) with respect to the burst period. However, this effect is likely mainly caused by an observational bias, that is, due to the fading of stellar populations with time low-mass clusters are more difficult to detect at old and intermediate age than at young ages (see also de Grijs et al. 2002).
In order to consider only cluster candidates with good photometry, de Grijs et al. (2001) selected clusters brighter than $`V=`$ 22.5 (corresponding to the 100% completeness limit). The equivalent upper limit in absolute visual magnitude is $`M_V=`$5.3. The temporal evolution of the corresponding lower mass limit is shown as the solid line in the two panels of Fig. 8. This fading line is based on the mass-to-light ratio predicted by the BC96 model for a metallicity of Z$``$ 0.02 Z. We can afford to use a model with such a low metallicity, and therefore a low mass-to-light ratio, since very metal-poor clusters are detected at any age up to the burst epoch. A comparison of the fading line with the distribution of points in top and bottom panels of Fig. 8 shows that the deficit in low-mass clusters at large age cannot be accounted for solely by the selection effect, especially in B1.

Systems of young clusters are well known to show a cluster initial mass function scaling as $`\mathrm{d}Nm^\alpha \mathrm{d}m`$, where $`m`$ is the cluster initial mass and $`\mathrm{d}N`$ is the number of clusters with mass between $`m`$ and $`m+\mathrm{d}m`$. Observations show that the slope $`\alpha `$ is $`2`$ (Zhang & Fall 1999; Whitmore et al. 1999; Bik et al. 2002; see also Sect. 4.6). At an intermediate age of 5 Gyr (log $`t9.7`$), the mass detection threshold is $`3\times 10^4\mathrm{M}_{}`$ (log $`m=4.5`$). If the cluster initial mass function had remained unaffected (above the fading line), the ratio between the number of clusters with masses $`3\times 10^4\mathrm{M}_{}\mathrm{m}3\times 10^5\mathrm{M}_{}(4.5\mathrm{log}m/\mathrm{M}_{}5.5)`$ and the number of clusters with masses $`3\times 10^5\mathrm{M}_{}m10^6\mathrm{M}_{}(5.5\mathrm{log}m/\mathrm{M}_{}6)`$ would be

$$\frac{_{0.03}^{0.3}m^2dm}{_{0.3}^1m^2dm}13.$$
(11)
To make use of the solar metallicity model raises the fading line by $`\mathrm{log}m0.3`$ at most in the $`\mathrm{log}m`$ vs $`\mathrm{log}t`$ plot through an increase of the metallicity dependent mass-to-light ratio. This therefore leads to a decrease in the number of detectable clusters whose mass is lower than $`3\times 10^5\mathrm{M}_{}`$. Even in that case however, the proportion between the number of high-mass clusters ($`5.5\mathrm{log}m/\mathrm{M}_{}6`$) and the number of lower-mass clusters (i.e. clusters whose mass is larger than the solar metallicity fading line and lower than $`3\times 10^5\mathrm{M}_{})`$ is still $`6.5`$. There is an obvious discrepancy between these theoretical ratios, whatever the metallicity of the fading line, and the observed ratios as the majority of the clusters in the pre-burst phase have masses greater than $`3\times 10^5\mathrm{M}_{}`$. We emphasize here that the mass error bars do not affect this outcome significantly. For instance, let us consider the old and high-mass cluster, B1-36 ($`\mathrm{log}t=10.`$, $`\mathrm{log}m/\mathrm{M}_{}=5.91`$; see Table 1). The rather large mass error bar suggest that its mass may be as low as $`\mathrm{log}m/\mathrm{M}_{}=4.86`$, so that it may contribute to the low-mass clusters of the pre-burst phase. However, its age and its mass errors are strongly correlated through the age-dependent mass-to-light ratio, that is, the lower mass limit is coupled with the lower age limit ($`\mathrm{log}t=8.6`$) and not with the best age estimate (see Sect. 4.2). As a result, the high-mass clusters whose minimal estimate is lower than $`\mathrm{log}m/\mathrm{M}_{}=5.5`$ are not expected to contribute to the low-mass clusters class of the pre-burst phase.

The depletion in low-mass clusters at intermediate age shown by Fig. 8 is a real effect and illustrates the dynamical evolution of the M82 B cluster system. The disruptive processes working on time-scales on the order of a Hubble time include the tidal interactions between the clusters and the gravitational field of the parent galaxy and the evaporation through two-body relaxation within clusters, the latter being enhanced by the former. They affect mostly the low-mass clusters (e.g. Gnedin & Ostriker 1997, de Grijs et al. 2002). The detection of such evolutionary effects is not unexpected because of the dense environment achieved in this fossil starburst of the low-mass galaxy M82. A similar cluster disruption time-scale was derived in the dense centre of the massive spiral M51 (Boutloukos & Lamers 2002). Furthermore, de Grijs et al. (2002) show that the mean densities in M82 B and in the centre of M51 are remarkably similar, within an order magnitude. This therefore suggests that cluster disruption time-scales are primarily driven by the density of the environment in which the clusters are embedded, irrespective of the mass of the host galaxy. Such a result is not unexpected. The cluster stars which venture to large enough distance from the cluster centre are stripped off due to the tidal influence of the host galaxy, leading to a cluster which is tidally truncated at a finite tidal radius $`r_t`$ given by:

$$r_t=R\left(\frac{m}{3M_{\mathrm{gal}}}\right)^{1/3},$$
(12)
where R is the galactocentric distance of the cluster and $`M_{\mathrm{gal}}`$ is the mass of the galaxy within $`R`$. Equation (12) also implies that the density of a GC, $`\rho _{\mathrm{cl}}`$, cannot be less than three times the average density of the host galaxy inside $`R`$, $`\rho _{\mathrm{gal}}`$, that is $`\rho _{\mathrm{cl}}3\rho _{\mathrm{gal}}`$. If its orbital motion drives the cluster towards higher density regions of the host galaxy, the corresponding shrinkage of the cluster tidal radius causes a runaway tidal stripping of the outer cluster stars, thus decreasing the cluster mass and making it even more sensitive to the galactic tidal field. Ultimately, this may lead to the cluster’s disruption.

Even in the Magellanic Clouds (MCs), which constitute environments dynamically and tidally more gentle than M82 B and M51, do low-mass clusters dissolve with time through normal internal process. Due to the smaller ambient density, the time-scale for depleting low-mass clusters is not that short in the MCs however: clusters with masses between $`10^4`$ and $`10^5M_{}`$ survive for a few Gyr (Elson & Fall 1985b, 1988, Mackey & Gilmore 2002a,b), which is not the case in M82 B. In the next section, we show that the depletion of low-mass clusters in M82 B is already at work at the burst time.

### 4.6 Mass functions

Figure 9 displays the cluster mass functions for B1 and B2 combined. We have considered the clusters with ages 8.6 $`<`$ log $`t/\mathrm{yr}`$ $`<`$ 9.1. To include older clusters in our sample would artificially bias the cluster mass function towards higher masses because of the dynamical destruction with time of the low-mass clusters. De Grijs et al. (2002) derive a characteristic cluster disruption time-scale for the M82 B cluster system of $`\mathrm{log}(t_{\mathrm{dis}}/\mathrm{yr})7.5+0.62\times \mathrm{log}(m/10^4\mathrm{M}_{}`$). Obviously, the dynamical destruction of star clusters in M82 B takes place on very short time-scales. In fact, the disruption time-scale in M82 B is the shortest known in any disc (region of a) galaxy. Considering the clusters formed at $`\mathrm{log}t/\mathrm{yr}=8.6`$, the relation between the disruption time-scale and the cluster mass mentioned above shows that significant disruption must have already affected clusters with masses below 5.9 $`\times `$ 105 M. For the clusters formed at $`\mathrm{log}t/\mathrm{yr}`$=9.1 (i.e. almost at the onset of the burst), the same relation indicates that the upper mass of the clusters significantly affected by disruption is even greater, i.e. 3.8 $`\times `$ 106 M. Thus, the dynamical destruction processes have strongly altered the whole range of the initial mass function of the M82 B clusters with ages 8.6 $`<`$ log $`t/\mathrm{yr}`$ $`<`$ 9.1. As a result, the current mass spectrum of these clusters (see Fig. 9) no longer represents the formation conditions of the cluster system.

In galaxies hosting systems of young clusters for which deep HST observations are available, the cluster luminosity functions do not show strong evidence for any turnover. In fact, their shape is consistent, down to the completeness threshold, with a decreasing power-law of the form $`\varphi (L)\mathrm{d}LL^\alpha \mathrm{d}L`$, where $`\varphi (L)\mathrm{d}L`$ is the number of young star clusters with luminosities between $`L`$ and $`L+\mathrm{d}L`$ and $`\alpha `$ is the slope of the cluster luminosity function. Since in the age range we consider, i.e. $`8.6<\mathrm{log}t/\mathrm{yr}<9.1`$, the mass-to-light ratio does not vary significantly with time and metallicity ($`m/L_V`$ $``$ 1), we can directly compare the mass function shown in Fig. 9 with the luminosity functions derived for other young star cluster systems in merging galaxies or merger remnants (e.g. NGC 4038/4039, $`\varphi (L)L^{1.8}`$, Whitmore & Schweizer 1995; NGC 3256, $`\varphi (L)L^{1.8}`$, Zepf et al. 1999; NGC 3921, $`\varphi (L)L^{2.1}`$, Schweizer et al. 1996; NGC 7252, $`\varphi (L)L^{1.8}`$, Miller et al. 1997; NGC 3610, $`\varphi (L)L^{1.9}`$, Whitmore et al. 2002). We have estimated the best-fitting power-law exponent with the Levenberg-Marquardt method (Press et al. 1997). The mass function in Fig. 9 follows $`\mathrm{d}Nm^{1.4\pm 0.2}\mathrm{d}m`$ (see also de Grijs et al. 2001). Our slope is therefore shallower than what has been found for many other young star cluster systems and shows better agreement with the luminosity function of the young star cluster population in the Large Magellanic Cloud ($`\alpha =1.5\pm 0.2`$, Elson & Fall 1985a). The most likely cause for such a rather shallow slope is the short disruption time-scale of the M82 B cluster system. The clusters first affected by the disruption are the lowest-mass ones, even though the most massive clusters formed at the time of the burst have already gone through significant disruption processes too. Due to the shorter disruption time-scale, any initial mass function will get biased towards higher masses as time goes on and an initial steep slope will turn into a shallower one. Moreover, the mass distribution of the clusters formed during the burst shows a turnover at $`\mathrm{log}m/M_{}5.3`$ (de Grijs et al. 2002), not caused by selection effects (see the location of the fading line in Fig. 8). This also contributes to the difference in slope.

## 5 Conclusions

We have estimated the age, the metallicity and the mass of 87 clusters located in M82 B, a fossil starburst site in the irregular galaxy M82. Since this galaxy is located outside the Local Group, our analysis is based on spectral synthesis models. We have compared the HST $`BVIJ`$ photometry obtained by de Grijs et al.(2001) with the colours of simple stellar populations (Bruzual & Charlot 1996). Because the isochrones and the iso-metallicity tracks are not parallel in a $`(VI)_0`$ vs $`(VJ)_0`$ plot, we have been able to lift the age-metallicity degeneracy.
Our results have confirmed the peak in the cluster formation rate detected by de Grijs et al. (2001). From 1.5 to 0.5 Gyr ago, as M82 was getting closer to its large spiral companion M81, M82 B went through a period of enhanced cluster formation, most probably induced by the increasing tidal interactions. We have emphasized its chemical evolution showing that M82 B underwent a chemical enrichment phase up to super-solar metallicities (i.e. Z$`>Z_{}`$) about 1 Gyr ago. The stellar ejecta from older stellar populations have most probably enriched the interstellar medium out of which the 1 Gyr old clusters formed. At almost the same time however, a subpopulation of very metal-poor clusters formed in B2, i.e. the part of M82 B nearest the galactic centre. Their formation may have been stimulated by infall, onto B2, of “fresh” circumgalactic gas which had escaped the chemical enrichment of the interstellar medium while orbiting M82. The clusters in B1, the eastern part of M82 B, may have been self-enriched at the time of their formation since they show a mass-metallicity correlation in the sense expected by simple self-enrichment models applied to pressure-truncated clouds. In B2 however, a clear-cut relation does not show up. Due to the infall of circumgalactic gas, the pressure may not have been stable enough to allow the formation of pressure-bound clouds. The eastern (B1) and western (B2) parts of M82 B are therefore different with respect to both the mass-metallicity diagrams and the formation of very metal-poor clusters at the time of the burst. We have also highlighted the dynamical destruction of the low-mass clusters over time-scales of the order of several $`10^7`$yr. Finally, we have shown that the mass function of the clusters in a given age range obeys a power-law as observed for many cluster systems formed in merging galaxies. However, the slope of the mass function of cluster systems associated with mergers is significantly steeper ($`\alpha 2`$) than the slope of the $``$ 1 Gyr old clusters in M82 B ($`\alpha 1.4`$). This is most likely due to the disruption processes at work in the M82 B region. Since the low-mass clusters are more quickly destroyed, the mass function will get biased towards higher masses and a steep initial mass function will turn into a shallower one.

## Acknowledgments

G.P. warmly thanks the Cambridge Institute of Astronomy for its hospitality where much of this work was carried out. This research was supported by the European Community under grant HPMT-CT-2000-00132. G.P. is grateful to P. Magain (Liège Institute) for fruitful advices when tackling the issue of the error propagation.