text:
# Chiral 1/𝑀² corrections to 𝐵^(∗)→𝐷^(∗) at Zero Recoil in Quenched Chiral Perturbation Theory

## I Introduction

The Cabbibo-Kobayashi-Maskawa (CKM) matrix describes the flavor mixing among the quarks, its elements are fundamental input parameters for the standard model. Their precise knowledge is not only crucial to determine the standard model but also to shed light on the origin of $`CP`$ violation. The matrix element that parameterizes the amount of mixing between the $`b`$ and $`c`$ quarks, $`V_{cb}`$, can be extracted from the exclusive semileptonic $`B`$ meson decays $`BDl\nu `$ and $`BD^{}l\nu `$, where $`l=e,\mu `$. Heavy quark effective theory (HQET) (for a recent review, see Manohar and Wise (2000)), which is exact in the limit of infinite masses $`M`$ for the heavy quarks, predicts the width of the process $`BD^{}l\nu `$ as

$$\frac{d\mathrm{\Gamma }}{d\omega }(BD^{})=\frac{G_F^2|V_{cb}|^2}{48\pi ^3}𝒦(\omega )_{BD^{}}(\omega )^2,$$
(1)
where $`\omega =vv^{}`$ is the scalar product of the 4-velocities $`v`$ and $`v^{}`$ of the $`D^{}`$ and $`B`$ mesons, respectively. $`𝒦(\omega )`$ is a known kinematical factor and $`(\omega )`$ is a form factor whose value at the kinematical point $`\omega =1`$ is $`(1)=1`$ in the $`M\mathrm{}`$ limit. There are, however, perturbative and non-perturbative corrections to $`(1)`$,

$$_{BD^{}}(1)=\eta _A+\delta _{1/M^2}+\mathrm{},$$
(2)
where the parameter $`\eta _A0.96`$ is a QCD radiative correction known to two-loop order Czarnecki (1996) and $`\delta _{1/M^2}`$ are non-perturbative corrections of $`𝒪(1/M^2)`$ to the infinite mass limit of HQET. Note that, according to Luke’s theorem Luke (1990) there are no $`𝒪(1/M)`$ corrections at zero-recoil. One chooses the zero-recoil point because, for $`\omega =1`$, $`_{BD^{}}`$ can be expressed in terms of a single form-factor $`h_{A_1}`$ given by

$$\frac{D^{}(v,ϵ^{})|\overline{c}\gamma ^\mu \gamma _5b|B(v)}{\sqrt{m_Bm_D^{}}}=2ih_{A_1}(1)ϵ^\mu .$$
(3)
This is in contrast to the general case $`\omega >1`$ for which $`_{BD^{}}(\omega )`$ is a linear combination of several different form factors of $`BD^{}l\nu `$ mediated by vector and axial vector currents.

Several experiments, most recently by CLEO Briere et al. (2002), have determined the product $`(_{BD^{}}(1)|V_{cb}|)^2`$ by measuring $`d\mathrm{\Gamma }_{BD^{}}/d\omega `$ and extrapolating it to the zero-recoil point. The mixing parameter $`|V_{cb}|`$ can then be extracted once the value $`_{BD^{}}(1)`$, that encodes the strong interaction physics, has been evaluated. The uncertainty in $`|V_{cb}|`$ is therefore determined by the experimental errors and by theoretical uncertainties in the determination of $`_{BD^{}}(1)`$. Presently, the theoretical uncertainties dominate.111 Similarly, one can use the decay $`BDl\nu `$ to extract $`(_{BD}(1)|V_{cb}|)^2`$ from the measured $`d\mathrm{\Gamma }_{BD}/d\omega `$. However, $`d\mathrm{\Gamma }_{BD}/d\omega `$ is more heavily suppressed by phase-space near $`\omega =1`$ than $`d\mathrm{\Gamma }_{BD^{}}/d\omega `$. In addition, the $`BD`$ channel is experimentally more challenging. Thus the extraction of $`|V_{cb}|`$ from this channel is less precise but serves as a consistency check.

A model-independent way of calculating $`(1)`$ is provided by numerical lattice QCD simulations. In this method one implements field theory non-perturbatively using the Feynman path integral approach. The fermion determinant that arises from the path integral is very costly to calculate; it is often set to one in an approximation called quenched QCD (QQCD). This corresponds to dropping the contribution from virtual quark loops which are made of “sea” quarks that have propagators not connected to the inserted external operators. “Valence” quarks, those that are connected to the inserted operators, however, are kept. Recently, such calculations have been performed Simone et al. (2000); El-Khadra et al. (2001); Hashimoto et al. (2001); Kronfeld et al. (2002) for the decays $`BD^{()}l\nu `$. Several systematic uncertainties, such as from statistics and lattice space dependence, contribute to the error of these calculations. Another contribution to the uncertainties comes from the chiral extrapolation of the light quark mass. Since lattice QCD simulations are limited by the available computing power they presently cannot be performed with the physical masses of the light quarks. Therefore one needs to extrapolate from the heavier masses used on the lattice (of order the strange quark mass) down to the physical light quark masses. This extrapolation can be done by matching QQCD to quenched chiral perturbation theory (Q$`\chi `$PT) and calculating the non-analytic corrections $`\delta _{1/M^2}`$ in Eq. (2) in Q$`\chi `$PT. The formally dominant contributions to these corrections come from the hyperfine mass splitting between the heavy pseudoscalar and vector mesons that stems from the inclusion of heavy quark symmetry breaking operators of $`𝒪(1/M)`$ in the Lagrangian.

In QCD, the corrections due to $`D`$ meson hyperfine splitting have been calculated in chiral perturbation theory ($`\chi `$PT) by Randall and Wise Randall and Wise (1993). A more complete treatment, involving additional corrections due to $`B`$ meson hyperfine splitting, $`𝒪(1/M)`$ axial coupling corrections, and $`𝒪(1/M)`$ corrections to the current, has been given in Boyd and Grinstein (1995a). Recently, the $`D`$ meson hyperfine splitting corrections have also been determined in partially quenched chiral perturbation theory (PQ$`\chi `$PT) Savage (2002) for partially quenched QCD (PQQCD). In contrast to QQCD, PQQCD does not drop the contributions from sea quark loops but gives different (and separately varied) masses to the sea and valence quarks. Usually, sea quarks are heavier than valence quarks, so that the fermion determinant — although no longer equal to one — is much less costly to calculate than in ordinary QCD.

In this paper we calculate the $`𝒪(1/M^2)`$ corrections in Eq. (2) due to $`D`$ and $`B`$ meson hyperfine splitting in quenched chiral perturbation theory (Q$`\chi `$PT). As we will show, these corrections are — upon expanding in powers of the hyperfine splitting $`\mathrm{\Delta }`$ — of order $`\mathrm{\Lambda }_{\text{QCD}}^{}{}_{}{}^{3n/2}/(M^nm_q^{n/2})`$ for $`n2`$ and formally larger than those coming from the inclusion of $`𝒪(1/M)`$ heavy quark symmetry breaking operators in the Lagrangian and current which are suppressed by powers of $`\mathrm{\Lambda }_{\text{QCD}}/M`$. This argument is similar to the one that applies to $`\chi `$PT Manohar and Wise (2000).

Our Q$`\chi `$PT calculation can be used to extrapolate lattice results Hashimoto et al. (2001) that use the quenched approximation down to the physical light quark masses. So far, this extrapolation has been based upon the $`\chi `$PT calculation Randall and Wise (1993). Using Q$`\chi `$PT should therefore give a better estimate of the uncertainties related to the chiral extrapolation.

A central role in the lattice calculation of $`BD^{}`$ Hashimoto et al. (2001); Kronfeld et al. (2002) is played by the double-ratios of matrix elements

$$_+=\frac{D|\overline{c}\gamma ^0b|BB|\overline{b}\gamma ^0c|D}{D|\overline{c}\gamma ^0c|DB|\overline{b}\gamma ^0b|B},$$
(4)
$$_1=\frac{D^{}|\overline{c}\gamma ^0b|B^{}B^{}|\overline{b}\gamma ^0c|D^{}}{D^{}|\overline{c}\gamma ^0c|D^{}B^{}|\overline{b}\gamma ^0b|B^{}},$$
(5)
and

$$_{A_1}=\frac{D^{}|\overline{c}\gamma ^j\gamma _5b|BB^{}|\overline{b}\gamma ^j\gamma _5c|D}{D^{}|\overline{c}\gamma ^j\gamma _5c|DB^{}|\overline{b}\gamma ^j\gamma _5b|B}.$$
(6)
Since the numerator and denominator are so similar, statistical fluctuations are highly correlated and cancel in the ratios to a large degree. The $`𝒪(1/M^2)`$ correction to the double ratios can therefore be calculated fairly accurately and used to derive the $`𝒪(1/M^2)`$ correction to the matrix elements themselves. For this reason, we also calculate $`𝒪(1/M^2)`$ corrections to the decay $`B^{}D^{}`$ in addition to the experimentally accessible decays $`BD`$ and $`BD^{}`$, and thus the corrections to $`_+`$, $`_1`$, and $`_{A_1}`$.

## II Quenched Chiral Perturbation Theory

In a world where all the quark masses are large compared to $`\mathrm{\Lambda }_{\text{QCD}}`$ internal quark loops are suppressed and the results from QQCD are close to those from QCD. In the real world, however, light quarks are light ($`\mathrm{\Lambda }_{\text{QCD}}`$) and contributions from internal quark loops are substantial. One can nevertheless study the low-energy behavior of QQCD by its effective low energy theory, Q$`\chi `$PT Morel (1987); Sharpe (1992); Bernard and Golterman (1992a, b); Golterman (1994).

We consider a theory constructed from three light valence-quarks, $`u`$, $`d`$, $`s`$, and three light bosonic quarks, $`\stackrel{~}{u}`$, $`\stackrel{~}{d}`$, $`\stackrel{~}{s}`$ governed by the Lagrangian of QQCD

$$=\underset{a=u,d,s}{}\overline{q}^a(iD/m_q)_a^aq_a+\underset{\stackrel{~}{a}=\stackrel{~}{u},\stackrel{~}{d},\stackrel{~}{s}}{}\overline{\stackrel{~}{q}}^{\stackrel{~}{a}}(iD/m_{\stackrel{~}{q}})_{\stackrel{~}{a}}^{\stackrel{~}{a}}\stackrel{~}{q}_{\stackrel{~}{a}}=\underset{j=u,d,s,\stackrel{~}{u},\stackrel{~}{d},\stackrel{~}{s}}{}\overline{Q}^j(iD/m_Q)_j^jQ_j.$$
(7)
Here both types of quarks have been accommodated in the six-component vector $`Q`$ with the three quarks $`q_a`$ in the upper three entries and the three bosonic ghost-quarks $`\stackrel{~}{q}_{\stackrel{~}{a}}`$ in the lower three entries. The graded equal-time commutation relation governing the valence- and ghost-quarks is

$$Q_i^\alpha (𝐱)Q_{j}^{\beta }{}_{}{}^{}(𝐲)(1)^{\eta _i\eta _j}Q_{j}^{\beta }{}_{}{}^{}(𝐲)Q_i^\alpha (𝐱)=\delta ^{\alpha \beta }\delta _{ij}\delta ^3(𝐱𝐲),$$
(8)
where $`\alpha `$ and $`\beta `$ are spin- and $`i`$ and $`j`$ are flavor-indices. The graded equal-time commutation relations for two $`Q`$’s and two $`Q^{}`$’s are analogous. $`\eta _k`$ is given by

$$\eta _k=\{\begin{array}{cc}1& \text{for }k=1,2,3\hfill \\ 0& \text{for }k=4,5,6\hfill \end{array}.$$
(9)
The quark mass matrix is given by $`m_Q=\text{diag}(m_u,m_d,m_s,m_u,m_d,m_s)`$, i.e., the fermionic and bosonic quarks have equal masses but different statistics. Therefore the contributions of fermionic and bosonic quarks in virtual quark loops cancel exactly. The Lagrangian in Eq. (7) exhibits a graded symmetry $`[SU(3|3)_LSU(3|3)_R]\times U(1)_V`$ that is assumed to be broken spontaneously to $`SU(3|3)_V\times U(1)_V`$. The dynamics of the emerging 36 pseudo-Goldstone mesons can be described at lowest order in the chiral expansion by the Lagrangian

$$=\frac{f^2}{8}\text{str}\left(^\mu \mathrm{\Sigma }^{}_\mu \mathrm{\Sigma }\right)+\lambda \text{str}\left(m_Q\mathrm{\Sigma }+m_Q^{}\mathrm{\Sigma }^{}\right)+\alpha ^\mu \mathrm{\Phi }_0_\mu \mathrm{\Phi }_0\mu _0^2\mathrm{\Phi }_0^2$$
(10)
where

$$\mathrm{\Sigma }=\mathrm{exp}\left(\frac{2i\mathrm{\Phi }}{f}\right)=\xi ^2$$
(11)
and

$$\mathrm{\Phi }=\left(\begin{array}{cc}\pi & \chi ^{}\\ \chi & \stackrel{~}{\pi }\end{array}\right).$$
(12)
Here the $`\pi `$, $`\stackrel{~}{\pi }`$, and $`\chi `$ are $`3\times 3`$ matrices of pseudo-Goldstone bosons with quantum numbers of $`\overline{q}q`$ pairs, pseudo-Goldstone bosons with quantum numbers of $`\overline{\stackrel{~}{q}}\stackrel{~}{q}`$ pairs, and pseudo-Goldstone fermions with quantum numbers of $`\overline{\stackrel{~}{q}}q`$ pairs, respectively,

$$\pi =\left(\begin{array}{ccc}\eta _u& \pi ^+& K^+\\ \pi ^{}& \eta _d& K^0\\ K^{}& \overline{K^0}& \eta _s\end{array}\right),\stackrel{~}{\pi }=\left(\begin{array}{ccc}\stackrel{~}{\eta }_u& \stackrel{~}{\pi }^+& \stackrel{~}{K}^+\\ \stackrel{~}{\pi }^{}& \stackrel{~}{\eta }_d& \stackrel{~}{K}^0\\ \stackrel{~}{K}^{}& \overline{\stackrel{~}{K^0}}& \stackrel{~}{\eta }_s\end{array}\right),\text{and}\chi =\left(\begin{array}{ccc}\chi _{\eta _u}& \chi _{\pi ^+}& \chi _{K^+}\\ \chi _\pi ^{}& \chi _{\eta _d}& \chi _{K^0}\\ \chi _K^{}& \chi _{\overline{K^0}}& \chi _{\eta _s}\end{array}\right).$$
(13)
The pion decay constant is fixed by Eq. (11) and $`f=132\text{ MeV}`$ in QCD.

The flavor-singlet field $`\mathrm{\Phi }_0`$ is defined as

$$\mathrm{\Phi }_0=\frac{1}{\sqrt{6}}\text{str}(\mathrm{\Phi })=\frac{1}{\sqrt{2}}(\eta ^{}\stackrel{~}{\eta }^{})$$
(14)
where str() denotes a supertrace over the flavor indices. $`\mathrm{\Phi }_0`$ is invariant under $`[SU(3|3)_LSU(3|3)_R]\times U(1)_V`$ and thus arbitrary functions of it can be included in the Lagrangian. To lowest order in the chiral expansion only the two operators included in Eq. (10) with parameters $`\alpha `$ and $`\mu _0`$ remain and are understood to be inserted perturbatively Bernard and Golterman (1992a). Notice that this singlet field $`\mathrm{\Phi }_0`$ is not heavy as in $`\chi `$PT and therefore cannot be integrated out. It introduces a new vertex, the so-called hairpin.

Upon expanding the Lagrangian in Eq. (10) one finds that the mesons with quark content $`u\overline{u}`$, $`d\overline{d}`$, and $`s\overline{s}`$, the only ones relevant for our calculation, have masses given by

$$m_{qq}^2=\frac{8\lambda m_q}{f^2}.$$
(15)
## III Inclusion of Heavy Quarks

The $`D`$-mesons with quantum numbers of $`c\overline{Q}`$ can be written as a six-component vector

$$D=(D_u,D_d,D_s,D_{\stackrel{~}{u}},D_{\stackrel{~}{d}},D_{\stackrel{~}{s}}).$$
(16)
Heavy quark symmetry is provided by combining creation and annihilation operators for the pseudoscalar and vector mesons, $`D`$ and $`D^{}`$ respectively, together into the field $`H^D`$

$`H^D`$ $`=`$ $`{\displaystyle \frac{1+v/}{2}}(D/^{}+i\gamma _5D),`$ (17)
$`\overline{H}^D=\gamma ^0H^D\gamma ^0`$ $`=`$ $`(D/_{}^{}{}_{}{}^{}+i\gamma _5D^{}){\displaystyle \frac{1+v/}{2}},`$ (18)
where $`v`$ denotes the velocity of a heavy meson. In HQET the momentum of a heavy quark is only changed by a small residual momentum of $`𝒪(\mathrm{\Lambda }_{\text{QCD}})`$. Hence, $`v`$ is not changed and $`H`$ is usually denoted by an index $`v`$ which we have dropped here to unclutter the formalism. In the heavy quark limit, the dynamics of the heavy mesons are described by the Lagrangian Booth (1995); Sharpe and Zhang (1996)

$$_D=i\text{tr}[\overline{H}_a^Dv_\mu (^\mu \delta _{ab}+iV_{ba}^\mu )H_b^D]+g\text{tr}(\overline{H}_a^DH_b^D\gamma _\nu \gamma _5A_{ba}^\nu )+\gamma \text{tr}(\overline{H}_a^DH_a^D\gamma _\mu \gamma _5)\text{str}A^\mu $$
(19)
where the traces tr() are over Dirac indices and supertraces str() over the flavor indices are implicit. The additional coupling term involving $`\mathrm{\Phi }_0\text{str}A^\mu `$ is a feature of Q$`\chi `$PT and not present in $`\chi `$PT. The light-meson fields are

$$A_\mu =\frac{i}{2}(\xi ^{}_\mu \xi \xi _\mu \xi ^{})=\frac{1}{f}_\mu \mathrm{\Phi }+𝒪(\mathrm{\Phi }^3)$$
(20)
and

$$V_\mu =\frac{i}{2}(\xi ^{}_\mu \xi +\xi _\mu \xi ^{})=\frac{i}{2f^2}[\mathrm{\Phi },_\mu \mathrm{\Phi }]+𝒪(\mathrm{\Phi }^4).$$
(21)
Expanding the Lagrangian $`_D`$ to lowest order in the meson fields leads to the (derivative) couplings $`DD^{}\varphi `$ and $`D^{}D^{}\varphi `$ whose coupling constants are equal as a consequence of heavy quark spin symmetry. The $`DD\varphi `$ coupling vanishes by parity.

An analogous formalism applies to the fields $`B`$ and $`B^{}`$ which are combined into $`H^B`$. Note that the axial coupling $`g`$ is the same for $`H^D`$ and $`H^B`$ mesons as dictated by heavy quark flavor symmetry.

We do not include terms of order $`m_q\sqrt{m_\pi }`$ in the Lagrangian as explicit chiral symmetry breaking effects are suppressed compared to the leading corrections. The presence of these terms is implied by the nonzero masses $`m_{qq}`$.

## IV Matrix elements of $`\overline{B}^{()}D^{()}l\overline{\nu }`$

The non-zero hadronic matrix elements for $`B^{()}D^{()}`$ can be defined in terms of the 16 independent form factors $`h_\pm `$, $`h_V`$, $`h_{A_{1,2,3}}`$, and $`h_{1\mathrm{}10}`$ as Falk and Neubert (1993); Manohar and Wise (2000)

$$\frac{D(v^{})|\overline{c}\gamma ^\mu b|B(v)}{\sqrt{m_Bm_D}}=h_+(\omega )(v+v^{})^\mu +h_{}(\omega )(vv^{})^\mu ,$$
(22)
$$\frac{D^{}(v^{},ϵ^{})|\overline{c}\gamma ^\mu b|B(v)}{\sqrt{m_Bm_D^{}}}=h_V(\omega )\epsilon ^{\mu \nu \alpha \beta }ϵ_\nu ^{}v_\alpha ^{}v_\beta ,$$
(23)
$$\frac{D^{}(v^{},ϵ^{})|\overline{c}\gamma ^\mu \gamma _5b|B(v)}{\sqrt{m_Bm_D^{}}}=ih_{A_1}(\omega )(\omega +1)ϵ^\mu +ih_{A_2}(\omega )(vϵ^{})v^\mu +ih_{A_3}(\omega )(vϵ^{})v^\mu ,$$
(24)
$`{\displaystyle \frac{D^{}(v^{},ϵ^{})|\overline{c}\gamma ^\mu b|B^{}(v,ϵ)}{\sqrt{m_B^{}m_D^{}}}}`$ $`=`$ $`(ϵ^{}ϵ)[h_1(\omega )(v+v^{})^\mu +h_2(\omega )(vv^{})^\mu ]+h_3(\omega )(ϵ^{}v)ϵ^\mu `$ (25)
$`+h_4(\omega )(ϵv^{})ϵ^\mu (ϵv^{})(ϵ^{}v)[h_5(\omega )v^\mu +h_6(\omega )v^\mu ],`$
and

$`{\displaystyle \frac{D^{}(v^{},ϵ^{})|\overline{c}\gamma ^\mu \gamma _5b|B^{}(v,ϵ)}{\sqrt{m_B^{}m_D^{}}}}`$ $`=`$ $`i\epsilon ^{\mu \alpha \kappa \delta }\{ϵ_\kappa ^{}ϵ_\delta [h_7(\omega )(v+v^{})^\mu +h_8(\omega )(vv^{})^\mu ]`$ (26)
$`+v^\alpha v^\beta [h_9(\omega )(ϵ^{}v)ϵ^\mu +h_{10}(\omega )(ϵv^{})ϵ^\mu ]\}.`$
Here, $`\omega =vv^{}`$ and $`v`$ \[$`ϵ`$\] and $`v^{}`$ \[$`ϵ^{}`$\] are the velocities \[polarization vectors\] of the initial state $`B^{()}`$ meson and final state $`D^{()}`$ meson, respectively. Note that we will not explicitly calculate matrix elements of $`B^{}D`$ as these can be easily related to the $`BD^{}`$ calculation by a hermitian conjugation of the matrix elements and an interchange of the $`c`$ and $`b`$ quarks, i.e., $`B^{()}D^{()}`$.

In the heavy quark limit the matrix elements in Eqs. (22)–(26) are reproduced by the operator

$$\overline{c}\gamma ^\mu (1\gamma _5)b\xi (\omega )\text{tr}[\overline{H}_v^{}^D\gamma ^\mu (1\gamma _5)H_v^B].$$
(27)
Here, $`\xi (\omega )`$ is the universal Isgur-Wise function Isgur and Wise (1989, 1990) with the normalization $`\xi (1)=1`$. To lowest order in the heavy quark expansion one finds

$$h_+(\omega )=h_V(\omega )=h_{A_1}(\omega )=h_{A_3}(\omega )=h_1(\omega )=h_3(\omega )=h_4(\omega )=h_7(\omega )=\xi (\omega )$$
(28)
and the remaining 8 form factors vanish.

The discussion of the $`\overline{B}^{()}D^{()}l\overline{\nu }`$ matrix elements is similar for different flavors of the light quark $`q`$ content of the $`B^{()}`$ and $`D^{()}`$ mesons; it applies equally to $`q=u`$, $`d`$, or $`s`$ as the theory splits into three similar copies of a one-flavor theory. In the limit of light quark $`SU(3)_V`$ flavor symmetry the matrix elements (and in particular the Isgur-Wise function) are therefore independent of the light quark flavor. However, in nature the masses of the $`u`$, $`d`$, and $`s`$ quarks are different and $`SU(3)_V`$ is not an exact symmetry. Therefore our results will include terms that depend upon $`m_q`$ via the meson masses $`m_{qq}`$ defined in Eq. (15).

## V $`1/M`$ Corrections

The lowest order heavy quark symmetry violating operator that can be included in the Lagrangian $`_D`$ in Eq. (19) is the dimension-three operator $`\frac{\lambda _{D_2}}{M_D}\text{tr}\left[\overline{H_a}^D\sigma ^{\mu \nu }H_a^D\sigma _{\mu \nu }\right]`$. It violates heavy-quark spin and flavor symmetries and comes from the QCD magnetic moment operator $`\overline{c}\sigma ^{\mu \nu }G_{\mu \nu }^AT^Ac`$, where $`G_{\mu \nu }^A`$ is the gluon field strength tensor and $`T^A`$ with $`A=1\mathrm{}8`$ are the eight color $`SU(3)`$ generators. This operator gives rise to a mass difference between the $`D`$ and $`D^{}`$ mesons of

$$\mathrm{\Delta }_D=m_D^{}m_D=8\frac{\lambda _{D_2}}{\overline{M}_D}.$$
(29)
This effect can be taken into account by modifying the $`D`$ and $`D^{}`$ propagators which become

$$\frac{i\delta _{ab}}{2(vk+3\mathrm{\Delta }_D/4+iϵ)}\text{and}\frac{i\delta _{ab}(g_{\mu \nu }v_\mu v_\nu )}{2(vk\mathrm{\Delta }_D/4+iϵ)},$$
(30)
respectively, so that in the rest frame, where $`v=(1,0,0,0)`$, an on-shell $`D`$ has residual energy of $`3\mathrm{\Delta }_D/4`$ and an on-shell $`D^{}`$ has residual energy of $`\mathrm{\Delta }_D/4`$. A similar effect due to the inclusion of a QCD magnetic moment operator for the $`b`$ quark applies to the $`B^{()}`$ mesons.

There are no corrections to the matrix elements for the semileptonic decays $`B^{()}D^{()}e\nu `$ of $`𝒪(1/M)`$ at zero-recoil according to Luke’s theorem Luke (1990). The leading corrections enter at $`𝒪(1/M^2)`$. In addition to tree-level contributions from the insertion of $`𝒪(1/M^2)`$ suppressed operators into the heavy quark Lagrangian or the current there are one-loop contributions from wavefunction renormalization and vertex correction. These one-loop diagrams have a non-analytic dependence on the meson mass $`m_{qq}`$ and depend on the subtraction point $`\mu `$. This dependence on $`\mu `$ is canceled by the tree-level contribution of the $`𝒪(1/M^2)`$ operators.

Because of the absence of disconnected quark loops in QQCD, which manifests itself as a cancellation between intermediate pseudo-Goldstone bosons and pseudo-Goldstone fermions in loops in Q$`\chi `$PT, the only loop diagrams that survive are those that contain a hairpin interaction or a $`\gamma `$ coupling.

The wavefunction renormalization contributions for the pseudoscalar and vector meson, $`Z_{D/B}`$ and $`Z_{D/B}^{}`$, respectively, come from the one-loop diagrams shown in Fig. (1) and have been calculated in Booth (1994, 1995); Sharpe and Zhang (1996). Including the $`\alpha `$ coupling we find for these diagrams

$`Z`$ $`=`$ $`1+{\displaystyle \frac{ig^2\mu _0^2}{f^2}}H_1(\mathrm{\Delta }){\displaystyle \frac{ig^2\alpha }{f^2}}H_2(\mathrm{\Delta })+{\displaystyle \frac{6i\gamma g}{f^2}}F_1(\mathrm{\Delta })`$ (31)
and

$`Z^{}`$ $`=`$ $`1+{\displaystyle \frac{ig^2\mu _0^2}{3f^2}}H_1(\mathrm{\Delta }){\displaystyle \frac{ig^2\alpha }{3f^2}}H_2(\mathrm{\Delta })+{\displaystyle \frac{2i\gamma g}{f^2}}F_1(\mathrm{\Delta })`$ (32)
$`+{\displaystyle \frac{2ig^2\mu _0^2}{3f^2}}H_1(0){\displaystyle \frac{2ig^2\alpha }{3f^2}}H_2(0)+{\displaystyle \frac{4i\gamma g}{f^2}}F_1(0).`$
The functions $`H_1`$, $`H_2`$, and $`F_1`$ come from loop integrals and are given in the Appendix. Note that in the heavy quark limit where $`\mathrm{\Delta }=0`$ one recovers $`Z=Z^{}`$, as required by heavy quark symmetry.

The vertex corrections come from one-loop diagrams. The non-vanishing contributions are shown in Fig. (2).

Combining the wavefunction renormalization and vertex corrections and including a local counterterm to cancel the dependence on the renormalization scale $`\mu `$ we find the following corrections for the form factors:

$`\delta h_+(1)`$ $`=`$ $`X_+(\mu )+{\displaystyle \frac{Z_B1}{2}}+{\displaystyle \frac{Z_D1}{2}}`$ (33)
$`{\displaystyle \frac{ig^2}{f^2}}\left[\mu _0^2H_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)\alpha H_8(\mathrm{\Delta }_B,\mathrm{\Delta }_D)\right]{\displaystyle \frac{6ig\gamma }{f^2}}G_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)`$
$``$ $`X_+(\mu )+{\displaystyle \frac{1}{(4\pi f)^2}}\left({\displaystyle \frac{g^2\mu _0^2}{3m^2}}\left[{\displaystyle \frac{g^2\alpha }{3}}2g\gamma \right]\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}\right)(\mathrm{\Delta }_B\mathrm{\Delta }_D)^2`$
$`+𝒪(\{\mathrm{\Delta }_B,\mathrm{\Delta }_D\}^3),`$
$`\delta h_{A_1}(1)`$ $`=`$ $`X_{A_1}(\mu )+{\displaystyle \frac{Z_B1}{2}}+{\displaystyle \frac{Z_D^{}1}{2}}`$ (34)
$`{\displaystyle \frac{ig^2}{3f^2}}\left[\mu _0^2H_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)\alpha H_8(\mathrm{\Delta }_B,\mathrm{\Delta }_D)\right]{\displaystyle \frac{2ig\gamma }{f^2}}G_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)`$
$`{\displaystyle \frac{2ig^2}{3f^2}}\left[\mu _0^2H_5(\mathrm{\Delta }_B,0)\alpha H_8(\mathrm{\Delta }_B,0)\right]{\displaystyle \frac{4ig\gamma }{f^2}}G_5(\mathrm{\Delta }_B,0)`$
$``$ $`X_{A_1}(\mu )+{\displaystyle \frac{1}{(4\pi f)^2}}\left({\displaystyle \frac{g^2\mu _0^2}{9m^2}}\left[{\displaystyle \frac{g^2\alpha }{9}}{\displaystyle \frac{2g\gamma }{3}}\right]\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}\right)\left(3\mathrm{\Delta }_B^2+\mathrm{\Delta }_D^2+2\mathrm{\Delta }_B\mathrm{\Delta }_D\right)`$
$`+𝒪(\{\mathrm{\Delta }_B,\mathrm{\Delta }_D\}^3),`$
and

$`\delta h_1(1)`$ $`=`$ $`X_1(\mu )+{\displaystyle \frac{Z_B^{}1}{2}}+{\displaystyle \frac{Z_D^{}1}{2}}`$ (35)
$`{\displaystyle \frac{ig^2}{3f^2}}\left[\mu _0^2H_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)\alpha H_8(\mathrm{\Delta }_B,\mathrm{\Delta }_D)\right]{\displaystyle \frac{2ig\gamma }{f^2}}G_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)`$
$`{\displaystyle \frac{2ig^2}{3f^2}}\left[\mu _0^2H_5(0,0)\alpha H_8(0,0)\right]{\displaystyle \frac{4ig\gamma }{f^2}}G_5(0,0)`$
$``$ $`X_1(\mu )+{\displaystyle \frac{1}{(4\pi f)^2}}\left({\displaystyle \frac{g^2\mu _0^2}{9m^2}}\left[{\displaystyle \frac{g^2\alpha }{9}}{\displaystyle \frac{2g\gamma }{3}}\right]\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}\right)\left(\mathrm{\Delta }_B\mathrm{\Delta }_D\right)^2`$
$`+𝒪(\{\mathrm{\Delta }_B,\mathrm{\Delta }_D\}^3),`$
which are defined by $`h_+(1)=1+\delta h_+(1)`$ and analog expressions for $`\delta h_{A_1}(1)`$ and $`\delta h_1(1)`$. The functions $`H_5`$, $`H_8`$, and $`G_5`$ come from loop-integrals that are listed in the Appendix and we have defined $`m=m_{qq}`$. The insertions of tree-level $`𝒪(1/M^2)`$ operators are represented by the functions $`X_+(\mu )`$, $`X_{A_1}(\mu )`$, and $`X_1(\mu )`$ which are independent of $`m`$ and exactly cancel the $`\mu `$ dependence of the logarithm. These functions can be extracted from lattice simulations by measuring the zero-recoil form factors for a varying mass of the light quark.

Experimentally, $`\mathrm{\Delta }_D142\text{MeV}`$ and $`\mathrm{\Delta }_B46\text{MeV}`$ so that the ratios $`\mathrm{\Delta }_D/m`$ and $`\mathrm{\Delta }_B/m`$, which enter the form factor corrections through the function $`R(\mathrm{\Delta }/m)`$ (defined in the Appendix), are $`𝒪(1)`$. On the lattice, however, one can vary all quark masses. Expanding first in powers of $`\mathrm{\Delta }`$ and then taking the chiral limit $`m0`$ one finds the formal limits given in Eqs. (33)–(35) where we have only kept the pieces non-analytic in $`m`$. This demonstrates that the terms linear in $`\mathrm{\Delta }_D`$ and $`\mathrm{\Delta }_B`$, although present in wavefunction renormalization and vertex corrections, cancel as required by Luke’s theorem Luke (1990). The leading order corrections are $`𝒪(\{\mathrm{\Delta }_B,\mathrm{\Delta }_D\}^2)`$.

As a consistency check one can restore heavy quark flavor symmetry by taking $`\mathrm{\Delta }_B=\mathrm{\Delta }_D`$. Since the $`𝒪(1/M^2)`$ corrections to $`h_+(1)`$ and $`h_1(1)`$ are proportional to $`(\mathrm{\Delta }_B\mathrm{\Delta }_D)^2`$ they disappear as they should since the charge associated with the operators $`\overline{c}\gamma _\mu c`$ and $`\overline{b}\gamma _\mu b`$ is conserved. This argument does not apply for the $`BD^{}`$ transition matrix element in the limit $`\mathrm{\Delta }_B=\mathrm{\Delta }_D`$ since there is no conserved axial charge associated with the operators $`\overline{c}\gamma _\mu \gamma _5c`$ and $`\overline{b}\gamma _\mu \gamma _5b`$.

In the chiral limit, the term proportional to $`\mu _0^2`$ has a $`1/m^2`$ singularity and dominates over the terms proportional to $`\alpha `$ and $`\gamma `$ that are only log-divergent. This is analogous to a term of the form $`(m_{qq}^2m_{jj}^2)/m_{qq}^2`$ found by Savage Savage (2002) for PQ$`\chi `$PT (here, $`m_{qq}`$ and $`m_{jj}`$ are valence and sea quark masses, respectively). In the limit $`m_{jj}m_{qq}`$ this term, however, vanishes as PQ$`\chi `$PT goes to $`\chi `$PT where the dominant term is $`\mathrm{log}m_{qq}`$. In Q$`\chi `$PT, on the other hand, the $`1/m_{qq}^2`$ pole persists, revealing the sickness of QQCD where the hairpin interactions give a completely different chiral behavior than in QCD.

The size of $`\mu _0`$ can be estimated from the $`\eta `$-$`\eta ^{}`$ mass splitting Sharpe (1992), large $`N_C`$ arguments Witten (1979); Veneziano (1979) ($`N_C`$ being the number of colors), or lattice calculations. These estimates imply $`\mu _0500900\text{MeV}`$; for the purpose of dimensional analysis we use $`\mu _0𝒪(\mathrm{\Lambda }_{\text{QCD}})`$. Taking $`g𝒪(1)`$ we therefore find that $`\delta h_+`$, $`\delta h_{A_1}`$, and $`\delta h_1`$ are of order $`\mathrm{\Delta }^n/m^n\mathrm{\Lambda }_{\text{QCD}}^{}{}_{}{}^{3n/2}/(M^nm_q^{n/2})`$ for $`n2`$ and thus larger than tree-level heavy quark symmetry breaking operators that are suppressed by $`\mathrm{\Lambda }_{\text{QCD}}/M`$.

To show the dependence of the zero-recoil form-factors on the mass of the light spectator quark it is necessary to know the numerical values of the parameters $`\mu _0`$, $`g`$, $`\alpha `$, and $`\gamma `$. In determining reasonable values for these couplings we follow the discussion by Sharpe and Zhang Sharpe and Zhang (1996). Assuming that $`g`$ is similar to the $`\chi `$PT value we use $`g^2=0.4`$. The hairpin coupling $`\alpha `$ is proportional to $`1/N_C`$, and thus assumed to be small; we use two values, $`\alpha =0`$ and $`\alpha =0.7`$. The coupling $`\gamma `$ is known to be suppressed by $`1/N_C`$ compared to $`g`$, the sign is undetermined. We take $`g\gamma g`$ (see Sharpe and Zhang (1996) and references therein).

With these parameters, the dependence of $`h_+(1)`$ and $`h_{A_1}(1)`$ on the mass of the light spectator quark $`m_q`$ is show in Figs. (3) and (4), respectively.

The graphs are plotted against $`m_q`$ in units of the strange quark mass $`m_s`$ with $`m_q/m_s=m^2/m_{\eta _s}^2`$ where $`m_{\eta _s}^2=2m_K^2`$. The behavior of $`h_+(1)`$ in Q$`\chi `$PT is dominated at small $`m`$ by the $`1/m^2`$ pole that is non-existent in $`\chi `$PT. Lattice calculations of $`h_+(1)`$ El-Khadra et al. (2001) show a small downward trend for decreasing $`m_q`$ down to the chiral limit that is similar to the downward trend seen from the $`\chi `$PT calculation (dashed line). The same behavior (down to $`m_q0.1m_s`$) can also be seen for Q$`\chi `$PT for a certain choice of parameters (e.g., $`\gamma `$ positiv). The case of $`h_{A_1}(1)`$ is different as there is a pole at $`m=\mathrm{\Delta }_D`$ which is close to the physical pion mass. Here, both $`D^{}`$ and $`\pi `$ can be on-shell and the decay $`BD^{}\pi `$ becomes kinematically allowed. Lattice calculations of $`h_{A_1}(1)`$ Hashimoto et al. (2001) for $`m_q=(0.6\mathrm{}1)m_s`$ show a small downward trend for decreasing $`m_q`$ similar to the downward trend seen from the $`\chi `$PT calculation \[dashed line in Fig. (4)\]. A similar trend down to $`m_q0.2m_s`$ can also be seen in the Q$`\chi `$PT calculation for a relatively large positiv value of $`\gamma `$.

Although the downward trend in the lattice data for the two cases seems significant as the statistical errors are highly correlated, the uncertainty is still relatively high (typically $`\pm 0.01`$) and the existing lattice data can be accommodated by a wide range of values for the parameters in the Q$`\chi `$PT Lagrangian.

Finally, we calculate the double ratios defined in Eqs. (4)–(6) using the results in Eqs. (33)–(35). We find

$$_+=1+2\delta h_+(1),$$
(36)
$$_1=1+2\delta h_1(1),$$
(37)
and

$`_{A_1}`$ $`=`$ $`1+\stackrel{~}{X}_{A_1}(\mu )`$ (38)
$`{\displaystyle \frac{ig^2}{3f^2}}\{\mu _0^2[H_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)+H_5(\mathrm{\Delta }_D,\mathrm{\Delta }_B)H_5(\mathrm{\Delta }_D,\mathrm{\Delta }_D)H_5(\mathrm{\Delta }_B,\mathrm{\Delta }_B)]`$
$`\alpha [H_8(\mathrm{\Delta }_B,\mathrm{\Delta }_D)+H_8(\mathrm{\Delta }_D,\mathrm{\Delta }_B)H_8(\mathrm{\Delta }_D,\mathrm{\Delta }_D)H_8(\mathrm{\Delta }_B,\mathrm{\Delta }_B)]\}`$
$`{\displaystyle \frac{2ig\gamma }{f^2}}\left[G_5(\mathrm{\Delta }_B,\mathrm{\Delta }_D)+G_5(\mathrm{\Delta }_D,\mathrm{\Delta }_B)G_5(\mathrm{\Delta }_D,\mathrm{\Delta }_D)G_5(\mathrm{\Delta }_B,\mathrm{\Delta }_B)\right]`$
$``$ $`1+\stackrel{~}{X}_{A_1}(\mu ){\displaystyle \frac{1}{(4\pi f)^2}}\left({\displaystyle \frac{2g^2\mu _0^2}{9m^2}}\left[{\displaystyle \frac{2g^2\alpha }{9}}{\displaystyle \frac{4g\gamma }{3}}\right]\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}\right)\left(\mathrm{\Delta }_B\mathrm{\Delta }_D\right)^2`$
$`+𝒪(\{\mathrm{\Delta }_B,\mathrm{\Delta }_D\}^3),`$
where $`\stackrel{~}{X}_{A_1}(\mu )`$ is the counter term associated with $`_{A_1}`$.

## VI Conclusions

Knowledge of the $`B^{()}D^{()}`$ form factors at the zero-recoil point is crucial to extract the value of $`V_{cb}`$ from experiment. In the limit that the heavy quarks are infinitely heavy HQET predicts that the form factors $`h_+`$, $`h_{A_1}`$, and $`h_1`$ are equal, $`h_+(1)=h_{A_1}(1)=h_1(1)=\xi (1)`$. The formally dominant correction due to breaking of heavy quark symmetry comes from the inclusion of a $`𝒪(1/M)`$ dimension-three operator in the Lagrangian that leads to hyperfine-splitting between the heavy pseudoscalar and vector mesons. These leading order corrections are $`𝒪(\{\mathrm{\Delta }_B,\mathrm{\Delta }_D\}^2`$ as required by Luke’s theorem.

Recent lattice simulations using the quenched approximation of QCD have made a big step forward in determining these zero-recoil form factors. Presently, however, the simulations use light quark masses that are much heavier than the physical ones and therefore rely on a chiral extrapolation down to the physical quark masses. In this paper we have calculated the dominant corrections to the form factors $`h_+`$, $`h_{A_1}`$, and $`h_1`$ in Q$`\chi `$PT and determined the non-analytic dependence on the light quark masses via the light meson masses $`m_{qq}`$. Using these results, instead of the $`\chi `$PT calculation, to extrapolate the QQCD lattice measurements of these form factors down to the physical pion mass should give a more reliable estimate of the errors associated with the chiral extrapolation.

We have also calculated the corrections to certain double ratios that are used in lattice QCD calculations of the decay $`BD^{}`$.

###### Acknowledgements.

I would like to thank Martin Savage for many helpful discussions. I am also grateful to him, Andreas Kronfeld, and Steve Sharpe for useful comments on the manuscript. This work is supported in part by the U.S. Department of Energy under Grant No. DE-FG03-97ER4014.

## Appendix A Integrals

We list the functions $`H_1`$, $`H_2`$, $`F_1`$, $`H_5`$, $`H_8`$, and $`G_5`$ even though some of them have appeared in the literature before Boyd and Grinstein (1995b, a). Here, $`m=m_{qq}`$ is the mass of the $`q\overline{q}`$ light meson in the loop where $`q=u`$, $`d`$, or $`s`$ is the light (spectator) quark content of the heavy mesons. We have used dimensional regularization with the $`\overline{\text{MS}}`$ scheme, where $`1/ϵ^{}1/ϵ\gamma _E+\mathrm{log}4\pi +1`$. At the end we set $`1/ϵ^{}=0`$. As a shorthand we have defined the function

$$R(x)=\sqrt{x^21}\mathrm{log}\left(\frac{x\sqrt{x^21+iϵ}}{x+\sqrt{x^21+iϵ}}\right),$$
(39)
which occurs frequently. We also need its derivative $`dR/dx`$ given by

$$R^{}(x)=\frac{x}{x^21}R(x)2.$$
(40)
For the calculation of the wavefunction renormalization contribution we need the derivatives of the loop integrals for the diagrams in Fig. (1):

$$H_1(\mathrm{\Delta })=\frac{i}{16\pi ^2}\left[\mathrm{log}\frac{m^2}{\mu ^2}\frac{1}{ϵ^{}}1R^{}\left(\frac{\mathrm{\Delta }}{m}\right)\right],$$
(41)
$`H_2(\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{i}{16\pi ^2}}[{\displaystyle \frac{16}{3}}\mathrm{\Delta }^2{\displaystyle \frac{10}{3}}m^2+2(m^2\mathrm{\Delta }^2)(\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}{\displaystyle \frac{1}{ϵ^{}}})+{\displaystyle \frac{4}{3}}\mathrm{\Delta }mR\left({\displaystyle \frac{\mathrm{\Delta }}{m}}\right)`$ (42)
$`+({\displaystyle \frac{2}{3}}\mathrm{\Delta }^2{\displaystyle \frac{5}{3}}m^2)R^{}\left({\displaystyle \frac{\mathrm{\Delta }}{m}}\right)],`$
and

$`F_1(\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{i}{16\pi ^2}}[{\displaystyle \frac{10}{3}}\mathrm{\Delta }^2{\displaystyle \frac{4}{3}}m^2+(m^22\mathrm{\Delta }^2)(\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}{\displaystyle \frac{1}{ϵ^{}}})+{\displaystyle \frac{4}{3}}\mathrm{\Delta }mR\left({\displaystyle \frac{\mathrm{\Delta }}{m}}\right)`$ (43)
$`+{\displaystyle \frac{2}{3}}(\mathrm{\Delta }^2m^2)R^{}\left({\displaystyle \frac{\mathrm{\Delta }}{m}}\right)].`$
For the loop-integrals of the vertex corrections one finds

$$H_5(\mathrm{\Delta },\stackrel{~}{\mathrm{\Delta }})=\frac{i}{16\pi ^2}\left(\mathrm{log}\frac{m^2}{\mu ^2}\frac{1}{ϵ^{}}1\frac{m}{\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }}}\left[R\left(\frac{\mathrm{\Delta }}{m}\right)R\left(\frac{\stackrel{~}{\mathrm{\Delta }}}{m}\right)\right]\right),$$
(44)
$`H_8(\mathrm{\Delta },\stackrel{~}{\mathrm{\Delta }})`$ $`=`$ $`{\displaystyle \frac{i}{16\pi ^2}}([2m^2{\displaystyle \frac{2}{3}}(\mathrm{\Delta }^2+\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{\mathrm{\Delta }}^2)](\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}{\displaystyle \frac{1}{ϵ^{}}})+{\displaystyle \frac{16}{9}}(\mathrm{\Delta }^2+\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{\mathrm{\Delta }}^2)`$ (45)
$`{\displaystyle \frac{10}{3}}m^2+{\displaystyle \frac{m(5m^22\stackrel{~}{\mathrm{\Delta }}^2)}{3(\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }})}}R\left({\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{m}}\right){\displaystyle \frac{m(5m^22\mathrm{\Delta }^2)}{3(\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }})}}R\left({\displaystyle \frac{\mathrm{\Delta }}{m}}\right)),`$
and

$`G_5(\mathrm{\Delta },\stackrel{~}{\mathrm{\Delta }})`$ $`=`$ $`{\displaystyle \frac{i}{16\pi ^2}}({\displaystyle \frac{10}{9}}(\mathrm{\Delta }^2+\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{\mathrm{\Delta }}^2){\displaystyle \frac{4}{3}}m^2+[m^2{\displaystyle \frac{2}{3}}(\mathrm{\Delta }^2+\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }}+\stackrel{~}{\mathrm{\Delta }}^2)](\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}}{\displaystyle \frac{1}{ϵ^{}}})`$ (46)
$`+{\displaystyle \frac{2m(\mathrm{\Delta }^2m^2)}{3(\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }})}}R\left({\displaystyle \frac{\mathrm{\Delta }}{m}}\right){\displaystyle \frac{2m(\stackrel{~}{\mathrm{\Delta }}^2m^2)}{3(\mathrm{\Delta }\stackrel{~}{\mathrm{\Delta }})}}R\left({\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{m}}\right)).`$


format: text

text:
# Clusters in the inner spiral arms of M51: the cluster IMF and the formation history Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555

## 1 Introduction

The interaction of galaxies triggers star formation as clearly seen in the young starbursts in interacting galaxies (for reviews see Kennicutt 1998; Schweizer 1998). One of the best examples is the Antennae system where the interaction is presently going on and star formation and cluster formation occurs on a very large scale in the region between the two merging nuclei (Whitmore et al. 1999).

The way the triggered star formation progresses through an interacting galaxy is best studied from a slightly older interacting system, where the close passage of the companion is over. In such a system one can hope to derive the age distribution of clusters as a function of location in the galaxy and thus measure the time sequence of the induced star and cluster formation as a function of location in the galaxy.

One of the best systems for this purpose is the interacting sytem of the Whirlpool galaxy (M51) at a distance of $`8.4\pm 0.6`$ Mpc (Feldmeier et al. 1997). The M51 system consists of a grand design spiral galaxy (NGC 5194) interacting with its dwarf companion (NGC 5195). The almost face-on orientation of M51 allows the observation of its structure in great detail with minimum obscuration by interstellar dust. The interaction of these two galaxies has been modelled by various authors, starting from the fundamental paper by Toomre & Toomre (1972). These authors derived an age of the closest passage to be $`2\times 10^8`$ years ago. Later refined models, e.g. by Hernquist (1990), and reviewed by Barnes (1998), have improved this time estimate. The closest approach is now believed to have occurred about 250 - 400 Myrs ago. The best model is found for a distance of the pericenter of 17 to 20 kpc, for a mass ratio of $`M/m2`$ and a relative orbit which crosses the plane of M51 under an angle of about 15 degrees. (For a discussion of the problems with this model and possible improvements, see Barnes 1998). Recently Salo & Laurikainen (2000) suggested that on the basis of $`N`$-body simulations that the M51 system had multiple passages, with the two last at 50 -100 Myrs and 400-500 Myrs ago.

Because of its relatively small distance from us, and the fact that we see the M51 system face one, the system is ideally suited for the study of the progression of cluster formation due to galaxy – galaxy interaction. For this reason we started a series of studies of different aspects of the M51 system based on *HST-WFPC2* observations in six broad band and two narrow band filters.

The nucleus of M51 was studied by Scuderi et al. (2002). They found that the core contains a starburst with an age of $`410\pm 140`$ Myrs and a total stellar mass of about $`2\times 10^7M_{\mathrm{}}`$ within the central 17 pc. This age agrees with the estimated time of closest passage of the companion, so the starburst in the core is most likely due to the interaction with the passing companion. M51 contains an unresolved nucleus with a diameter smaller than 2 pc and a luminosity of $`2\times 10^6L_{\mathrm{}}`$ (Scuderi et al. 2002).

The bulge, i.e. the reddish region with a size of $`11\times 16`$ arcsec2 = $`460\times 680`$ pc2, between the nucleus and the inner spiral arms, was studied by Scuderi et al. (2002) and by Lamers et al. (2002). The bulge is dominated by an old stellar population with an age in excess of 5 Gyrs. The $`HSTWFPC2`$ images of the bulge show the presence of dust lanes. The total amount of dust in the bulge is about $`2.3\times 10^3`$ $`M_{}`$ and the dust has about the same extinction law and approximately the same gas to dust ratio as our Galaxy (Lamers et al. 2002). This suggests a metallicity close to that in the solar neighbourhood. The $`HSTWFPC2`$ images of the bulge of M51 show clearly that the dust is concentrated in structures in the forms of spiral-like dust lanes and a bar that reaches all the way down into the core. Most intriguingly is the discovery of about 30 bright and mainly blue point-like sources in the bulge that are aligned more or less along the spiral-like dust lanes. Lamers et al. (2002) have shown that these are most likely very young massive stars $`20<M_{}<150M_{}`$ with little evidence of associated clusters. This mode of star formation is the result of the peculiar conditions, in particular the destruction of CO molecules, of the interstellar clouds in the bulge of M51.

In this study we focus on the clusters at a distance of about 1 to 3 kpc from the nucleus, i.e. near the inner spiral arms. The purpose of the paper is two-fold:
– (a) to determine the cluster intitial mass function, and
– (b) to determine the presence or absence of a starburst period that can be linked to triggering by the passage of the companion. This is not an easy task, because our data will show that the age distribution of the clusters is strongly affected by the disruption of clusters older than about 40 Myrs.

We study the clusters and their properties by identifying point-like sources in the $`HSTWFPC2`$ images and measuring their $`UBVRI`$ magnitudes to obtain their energy distributions. The magnitudes indicate that the sources are clusters instead of single stars. Their energy distributions are compared to cluster evolutionary synthesis models to determine the age, mass and $`E(BV)`$ of these clusters.

In Sect. 2 we describe the observations and the data reduction. In Sect. 3 the selection of cluster candidates is discussed. In Sect. 4 we describe the cluster evolutionary synthesis models and in Sect. 5 the fitting procedure for the derivation of the cluster parameters is explained. The mass versus age distribution of the clusters is derived in Sect. 6. In Sect. 7 the initial mass function of the clusters is derived from the sample of clusters younger than 10 Myrs. The cluster formation history is studied in Sect. 8. The summary and conclusions are given in Sect. 9.

## 2 Observations and reduction

M51 was observed with *HST-WFPC2* as part of the *HST Supernova INtensive Study (SINS)* program (Millard et al. 1999). For this study we use the images taken in the broad band filters F336W ($`U`$), F439W ($`B`$), F555W ($`V`$), F675W ($`R`$) and F814W ($`I`$) from the SINS program and in the narrow band filters F502N $`(\text{[OIII]})`$ and F656N (H$`\alpha `$) from the GO-program of H.C. Ford. The image in $`U`$ was taken on 1994 May 12, the $`BVRI`$ images were taken on Jan 15 1995 and the \[OIII\] and H$`\alpha `$ images on Jan 25 1995. The $`U`$ and $`B`$ images were split into three and two exposures of 400 $`s`$ and 700 $`s`$ respectively. The \[OIII\] and H$`\alpha `$ images are split into two exposures of 1200 $`s`$ and 500 $`s`$ (\[OIII\]), and 1400 $`s`$ and 400 $`s`$ (H$`\alpha `$). In the remaining bands one single exposure of 600 $`s`$ was taken. The data was processed through the PODPS (Post Observing Data Processing System) for bias removal, flat fielding and dark frame correction.

To remove the cosmic rays from the $`U`$, $`B`$ \[OIII\] and H$`\alpha `$ images, we used the STSDAS task *crrej* for combining the available exposures. For the $`VRI`$ images, where only one exposure is available, we used a procedure called “Cosmic Eraser”. This procedure combines the IRAF tasks *cosmicrays* and *imedit* to reject as carefully as possible the cosmic rays. The automatic detection of cosmic rays with the task *cosmicrays* is based upon two parameters, a detection threshold and a flux ratio. The first parameter enables the detection of all the pixels with a value larger than the average value of the surrounding pixels. The flux ratio is defined as the percentage of the average value of the four neighboring pixels (excluding the second brightest pixel) to the flux of the brightest pixel. This parameter allows a classification of the detected objects: cosmic ray or star. Training objects are used to determine the flux ratio carefully. These training objects are labeled by the user to be a cosmic ray or a star. With *imedit* the detected cosmic rays signals are replaced by an interpolation of a third order surface fit to the surrounding pixels.

After the correction for the cosmic rays, the images were corrected for bad pixels using the hot pixel list from the STScI $`WFPC2`$ website in combination with the task *warmpix*. Corrections for non-optimal charge transfer efficiency on the CCD’s of the $`WFPC2`$ camera were applied using the formulae by Whitmore & Heyer (1997).

With the task *daofind* from the DAOPHOT package Stetson (1987), we identified the point sources on the image. We performed aperture photometry on these sources, also with the DAOPHOT package. We used an aperture radius of 3 pixels. The sky background was calculated in an annulus with internal and external radius of 10 and 14 pixels respectively. We only selected the point sources with an uncertainty smaller than 0.2 in the magnitude. Photometric zeropoints were obtained from table 28.2 of the *HST Data Handbook* (Voit 1997), using the VEGAMAG photometric system (Holtzman et al. 1995).

The aperture correction was measured for a number of isolated, high S/N point sources on each $`WFPC2`$-chip. This output was adopted for all the other point sources on the chip. Following Holtzman et al. (1995) we have normalized the aperture correction to 1″(10 WF pixels). The aperture corrections we found are between -0.24 and -0.37 magnitudes. This is larger than the aperture corrections for stars $``$ -0.17 magnitudes (Holtzman et al. 1995), which means that the detected point sources are fairly well resolved.

We adopt a distance of $`d=8.4\pm 0.6`$ Mpc (Feldmeier et al. 1997), which corresponds to a distance modulus of 29.62. At this distance, 1 ″ corresponds to a linear distance of 40.7 pc, which means that an $`HSTWFC`$ pixel of 0.1″ corresponds to 4.1 pc.

## 3 Selection of the clusters

The $`U`$-band images were taken in 1994, 8 months earlier than the images in the other passbands. The orientation of $`HSTWFPC2`$ was not the same at the two epochs. In the 1995 $`BVRI`$-images the nucleus of M51 was in the center of the PC-image. The $`U`$-image was centered on SN1994I (that was 15” off from the nucleus) so that the nucleus was near the edge of the PC-image. The orientation of the $`WFPC2`$-images is shown in Fig. 1.

For our analysis only the part of $`WF4`$-chip of the $`U`$ image that overlap with the $`WF2`$-chip of the $`BVRI`$ images is used (see Fig. 1). To identify the various sources, we obtained a position of the sources in $`U`$ and transformed these positions to the coordinate frame of the $`BVRI`$ images. The orientation of the \[OIII\] and H$`\alpha `$ images is very similar to those of the $`BVRI`$ images. We transformed these narrow band images to the same orientation as the $`BVRI`$ images.

We applied our analysis of the cluster energy distributions to the sources that were detected in the images of at least three of the $`BVRI`$ broad-band filters with a magnitude uncertainty smaller than 0.2 in each band. To this purpose we compared the positions of all detected sources in the $`BVRI`$ images. If objects occur in at least three images within a position tolerance of 2 pixels they are selected for our study. This tolerance was chosen on the basis of tests with clearly identifiable sources that were observed in most bands. For the sources detected in at least three images we then checked for the presence of a source in the $`U`$ image, in the region where the $`U`$ image and the $`BVRI`$ images do overlap (the dark grey region in Fig. 1). If a source was not detected in a particular band, we adopted a conservative lower limit for the magnitude in that band (see below). This was not always possible for the $`U`$-magnitude, because many of the sources detected in the $`BVRI`$ bands are located in an area that was not covered by the $`U`$-image. As a consequence, we have no information on the $`U`$-magnitude for many of the sources.

The location of the sources in the $`WF2`$-chip is shown in Fig. 2. The overlay of the sources on the $`V`$ band image clearly shows that the clusters are concentrated in or near the spiral arms at a distance of about 1.5 kpc from the nucleus of M51. The number of sources, detected in the various bands, is listed in Table 1. The sample contains a total of 877 sources for which we have reliable photometry ($`\sigma <0.20`$ magn.) in at least three bands. 111The coordinates and the photometry of the objects are available in the electronic version of A& A

Fig. 3 shows the magnitude distribution of the sources in the various bands. All distributions show a slow increase in numbers towards fainter objects, a maximum and a steep decrease towards even fainter sources. The slow increase to fainter sources reflects their luminosity function. The steep decrease to high magnitudes is due to the detection limit. This detection limit is not a single value for each band, because of the variable background in the $`WF2`$ field due to the spiral arms. The maximum of the distributions are near $`U20.0`$, $`B22.0`$, $`V22.0`$, $`R21.5`$ and $`I21.0`$. These values will be adopted as conservative lower limits for the magnitudes of sources that were not detected in any given band.

All the objects fall in the visual magnitude range of $`16.5<V<23.6`$. For distance modulus of 29.62 this corresponds to an absolute magnitude range of $`12.1<M_V<6.0`$ in the case of zero extinction and $`12.6<M_V<6.5`$ in case of moderate extinction with $`E(BV)0.17`$. This means that the faintest objects could be either very bright stars or small clusters. The vast majority of the selected point sources are so bright that they must be clusters. (See Sect 5.6).

### 3.1 The H$`\alpha `$ and \[OIII\] magnitudes

The magnitudes in the narrow-band images of filters $`F656N`$ ($`m(H\alpha )`$) and $`F502N`$ ($`m[OIII]`$) can be compared with those of the wide-band filters $`F675W`$ ($`R`$) and $`F555W`$ ($`V`$) respectively to derive a measure of the equivalent width of the H$`\alpha `$ and \[OIII\] lines. When the magnitudes are expressed in the ST-system, the magnitude differences $`Rm(H\alpha )`$ and $`Vm[OIII]`$ are by definition equal to 0.0 if the spectrum of a source is a featureless, flat continuum (i.e. the spectrum of the source is a continuum $`F_\lambda =const`$ without a photospheric H$`\alpha `$ absorption, and there is no HII region around the star). The differences $`Rm(H\alpha )`$ and $`Vm[OIII]`$ are either positive or negative if the source spectrum has an emission or an absorption line, respectively, falling in the narrow band filter, independently of the interstellar extinction. With the Vega-system magnitudes one has to fold in both the non-zero spectral slope and the discrete features present in $`\alpha `$ Lyrae’s spectrum. Therefore, a featureless flat continuum would correspond to $`Rm(H\alpha )+0.08`$ and $`Vm[OIII]0.26`$ colors in the Vega-system.

Fig. 4a shows the $`m(H\alpha )`$ versus $`R`$ and Fig. 4b shows $`m[OIII]`$ versus $`V`$ for the subset of sources for which we could measure these magnitudes with an accuracy better than 0.2. The figure shows that many sources have H$`\alpha `$ emission ($`m(H\alpha )<R`$) but none (!) of the sources has \[OIII\] emission (note that $`m[OIII]`$ is slightly fainter than the $`V`$ magnitude as appropriate for complete absence of an emission line). The lack of \[OIII\] emission indicates that the far-UV flux, at $`\lambda <350`$ Å i.e. at energies higher than the second ionization potential of oxygen, is quite low. For main sequence stars with $`L>2.5\times 10^5L_{}`$ (i.e. with $`M>30M_{}`$ and $`T_{\mathrm{eff}}>\mathrm{33\; 000}`$ K), surrounded by an HII region with approximately solar abundances, we would expect $`Vm[OIII]>Rm(H\alpha )`$. We will show later, on the basis of the study of the energy distributions, that many of the objects are very young clusters with ages less than a 10 Myr. The lack of \[OIII\] emission shows that these clusters do not contain stars with $`T_{\mathrm{eff}}`$ above about 30 000 K, which corresponds to about 25 to 30 $`M_{}`$ (e.g. Chiosi & Maeder 1986). We can exclude the possibility that the lack of \[OIII\] emission be due to just a high metallicity, which could reduce the electron temperature in an HII region and weaken optical forbidden lines considerably, without requiring a lower effective temperature, or, equivalently, a low upper mass cutoff. This is because in many HII regions in M51 an in other metal-rich spiral galaxies not only are the \[OIII\] lines faint or absent but also the HeI lines are unusually faint relative to Balmer lines (Panagia 2000; Lenzuni & Panagia, in preparation). This result indicates that He is only partially ionized and, therefore, that the ionizing radiation field is indeed produced exclusively by stars with effective temperatures much lower than 33,000K . So we conclude that the upper mass limit for clusters in the inner spiral arms of M51 is about 25 to 30 $`M_{}`$.

## 4 Cluster evolution models

To determine the age, initial mass and $`E(BV)`$ of the clusters, we compare the observed energy distributions with energy distributions derived from theoretical cluster evolutionary synthesis models. We used two sets of models: the Starburst99 models for ages up to 1 Gyr, and the Frascati models for ages of 10 Myr to 5 Gyr. In this section the models are described.

### 4.1 Starburst99 cluster evolution models

To determine the age, initial mass and $`E(BV)`$ of the clusters, we compare the observed energy distributions with energy distributions derived from theoretical cluster evolutionary synthesis models. We used two sets of models: the Starburst99 models for ages up to 1 Gyr, and the Frascati models for ages of 10 Myr to 5 Gyr. In this section the models are described.

### 4.2 Starburst99 cluster evolution models

We compared the observed spectral energy distributions (SED) of our objects with the cluster evolutionary synthesis models of Leitherer et al. (1999), i.e. the Starburst99 models for instantaneous star formation. In these models, the stellar atmosphere models of Lejeune et al. (1997) are included. For stars with a strong mass loss the Schmutz et al. (1992) extended model atmospheres are used. The stellar evolution models of the Geneva group are included in these models. For our analysis we used the models of the Spectral Energy Distributions (SEDs).

In the cluster models the stars are assumed to be formed with a classical Salpeter IMF-slope of $`d\mathrm{log}(N)/d\mathrm{log}(M)=2.35`$. The lower cut-off mass is $`M_{\mathrm{low}}`$ = 1 M. We adopt two values for the upper cut-off mass: the standard value of $`M_{\mathrm{up}}`$ = 100 M, and the value of $`M_{\mathrm{up}}=30M_{}`$. This last value is adopted because it is suggested by the lack of \[OIII\] emission (see 3.1). The total initial mass of all models is $`10^6`$ $`M_{}`$. Leitherer et al. (1999) present models for 5 different metallicities, from $`Z=0.001`$ to $`Z=0.040=2\times Z_{\mathrm{}}`$. Observations of HII-regions have shown that the metallicity of the inner region of M51 is approximately $`2Z_{}`$ or slightly higher (e.g. Diaz et al. 1991; Hill et al. 1997). Therefore we adopt the models with $`Z=0.020`$ and 0.040. The almost solar metallicity agrees with the study of Lamers et al. (2002), who found that the extinction properties and the gas to dust-ratio in the bulge of M51 is very similar to that in the Milky Way. So we use four sets of models, denoted by the pair $`(Z,M_{\mathrm{up}})`$ of (0.02,100), (0.02,30), (0.04,100) and (0.04,30). For these combinations we calculated the SEDs of the cluster models, taking into account the nebular continuum emission.

To obtain the absolute magnitudes from the theoretical spectral energy distributions in the five broad-band filters and the two narrow-band filters we convolved the predicted spectra of the Starburst99 models with the $`WFPC2`$ filter profiles. The spectrum of Vega was also convolved with the filter functions in order to find the zero-points of the magnitudes in the VEGAMAG system. For each combination of $`Z`$ and $`M_{\mathrm{up}}`$ we calculated the SED of 195 models in the range of 0.1 Myr to 1 Gyr, with time steps increasing from 0.5 Myr for the youngest models to 10 Myr for the oldest models.

The adopted lower mass limit of 1 $`M_{}`$ for the cluster stars has consequences for the derived masses of the clusters. If the IMF with a slope of -2.35 has a lower mass limit of 0.2 $`M_{}`$, as for the Orion Nebula Cluster (Hillenbrand 1997), the derived masses of the clusters would be a factor 2.09 higher for an upper limit of 30 $`M_{}`$. If the lower limit is 0.6 $`M_{}`$, the mass of the clusters will be a factor 1.28 higher than those of the Starburst99 models for $`M_{up}=30`$ $`M_{}`$. We will take this effect into account in the determination of the cluster masses.

Some of the theoretical energy distributions in the $`WFPC2`$ wide-band filters used in this study are shown in Fig. 5, for clusters with an initial mass of $`10^6`$ $`M_{}`$ at the distance of M51. During the first 5 Myrs the UV flux remains high because most of the O-type stars have main sequence ages longer than 5 Myr. Between 5 and 10 Myr the O-type stars disappear. This results in a general decrease at all magnitudes, except in the $`I`$ band where the red supergiants start to dominate. At later ages the flux decreases at all wavelengths as stars end their lives. After 900 Myrs only stars of types late-B and later have survived.

Fig. 5 shows that at any wavelength in the observed range the flux decreases with time. This fading of the clusters indicates that the detection limit of the clusters will gradually shift to those of higher mass as the clusters age. This is because the flux of a cluster (for a given stellar initial mass function) at any age and wavelength is proportional to the initial number of stars. The more massive the initial cluster, the brighter the cluster and the longer the flux will remain above the detection limit as the cluster fades due to stellar evolution. For any given detection limit, we can calculate the lower mass limit of the observable cluster as a function of age.

Fig. 6 shows the relation between the initial mass of a cluster and its age when it fades below the detection limit, derived from the Starburst99 models. The $`R`$-magnitude ($`R`$) of a cluster of initial mass $`M_i`$ without extinction at the adopted distance of 8.4 Mpc of M51 ($`dm=29.62`$) is related to the absolute $`R`$-magnitude ($`M_R`$) predicted for the Starburst99 models of $`10^6`$ $`M_{}`$ by

$$R(t)=M_R(t)+29.622.5\times \mathrm{log}(M/10^6)$$
(1)
So, for a given detection limit $`R_{\mathrm{lim}}`$ the critical mass $`M_{\mathrm{lim}}`$ of a cluster that can be detected at age $`t`$ is given by

$$\mathrm{log}M_{\mathrm{lim}}(t)=6+0.4\times (M_R(t)+29.62R_{\mathrm{lim}})$$
(2)
We have adopted the $`R`$ magnitudes to calculate the observed lower limits of the distributions for two reasons: (a) all clusters were detected in the $`R`$ band and (b) the effect of extinction is small in the $`R`$ band. The resulting relation between $`M_{\mathrm{lim}}`$ and age is shown in Fig. 6 for detection limits of $`R_{\mathrm{lim}}`$= 22.0 and 23.0 magn. We will use these limits later, in Sect. 5, to explain the distribution of clusters in the mass versus age diagram.

Fig. 7 shows the variations of the colours of the Starburst99 models for solar metallicity and with an upper mass limit of 30 $`M_{}`$. The figure shows that the colours do not change monotonically, but that there are peaks and dips, especially in the age range of $`6.5<\mathrm{log}(t)<7.3`$. These are related to the phases of the appearance of red supergiants that also produced the dips in the limiting magnitude curves of Fig. 6.

### 4.3 The Frascati models

Since the Starburst99 models only cover the age range up to 1 Gyr, we also used synthetic cluster models for older ages from the Frascati group. These ”Frascati-models” were calculated by Romaniello (1998) from the evolutionary tracks of Brocato & Castellani (1993) and Cassisi et al. (1994) using the $`HSTWFPC2`$ magnitudes derived from the stellar atmosphere models by Kurucz (1993). These models are for instantaneous formation of a cluster of solar metallicity stars in the mass range of 0.6 to 25 $`M_{}`$, with a total initial mass of 50 000 $`M_{}`$, distributed according to Salpeter’s IMF. These models cover an age range of 10 to 5000 Myr. There are 48 models with timesteps increasing from 10 Myr for the youngest ones to 500 Myr for the oldest ones. We have increased the magnitudes of the Frascati models by -3.526 magn. in order to scale these models to a total mass of $`10^6`$ $`M_{}`$ in the mass range of 1.0 to 25 $`M_{}`$; i.e. $``$3.25 magn. correction for the conversion from $`\mathrm{5\; 10}^4`$ to $`\mathrm{1\; 10}^6`$ $`M_{}`$, and 0.279 magn. correction for the conversion from a lower stellar mass limit of 0.6 $`M_{}`$ to 1.0 $`M_{}`$. In this way the masses can be compared directly with those of the Starburst99 models.

The Starburst99 models are expected to more accurate for the younger clusters, because they are based on the evolutionary tracks of the Geneva-group which include massive stars. The Frascati models are expected to be more accurate for the old clusters, because they include stars with masses down to 0.6 $`M_{}`$.

### 4.4 The extinction curve

To determine the values of $`E(BV)`$ of the clusters we assume the mean galactic reddening law. Lamers et al. (2002) and Scuderi et al. (2002) have shown that the extinction law of respectively the bulge and the core of M51 agree very well with the galactic law. The values of $`R_i=A_i/E(BV)`$ for the $`HSTWFPC2`$ filters have been calculated by Romaniello (1998). They are respectively 4.97, 4.13, 3.11, 2.41 and 1.91 for the $`U,B,V,R`$ and $`I`$ bands.

In the fitting of the energy distributions we did not take into account the disruption of clusters. We assumed that the stars contribute to the total energy distribution of the clusters until the cluster is completely dissolved and is no longer detectable as a point source. This is a reasonable assumption because we measured the magnitudes in a circle of 40.7 pc radius (see Sect. 2).

## 5 Fitting the observed energy distributions to observed cluster models

### 5.1 The fitting procedure

We have fitted the energy distributions of the clusters with good photometry in four or more filters, viz. $`UBVRI`$, $`BVRI`$, $`UVRI`$ and $`UBVR`$ (see Table 1), with the energy distributions of the cluster evolutionary synthesis models, discussed above, using a three dimensional maximum likelihood method. The three fitting parameters are: the age of the cluster ($`\tau `$), the reddening ($`E(BV)`$) and the initial mass of the cluster ($`M`$). For the clusters detected in only three filters we reduce the parameter space and make a two dimensional maximum likelihood fit. For objects not detected in a band, we adopted the magnitude lower limits (described in Sect. 3) to check which fits were acceptable.222We will refer to this method in subsequent papers as the ”3(2)DEF- method”, i.e. the three (or two) dimensional energy distribution fitting”.

### 5.2 Three dimensional maximum likelihood method

This method was applied to sources that are detected in at least four bands. We fitted the observed energy distributions with those predicted for the Starburst99 cluster models and the Frascati models (see Sect. 4.2). For each model-age we have two fitting parameters $`M`$ and $`E(BV)`$. The cluster models are for an initial cluster mass of $`10^6`$ $`M_{}`$. For other masses the flux simply scales $`M/10^6M_{}`$. We have adopted an uncertainty in the model fluxes of 5 percent (0.05 magnitudes) in all bands.

To reduce the range in masses in the parameter space, we make an initial guess of the cluster mass based on the observed magnitudes, since the distance to M51 is known. By calculating the average difference in magnitude – weighted by the errors – between the theoretical energy distribution for $`10^6`$ $`M_{}`$ and the observed one, we have a good first approximation of the mass of the cluster:

$`M_{guess}(E(BV),\tau )`$
$`{\displaystyle \frac{_\lambda 1/\sigma _\lambda \{m_\lambda ^{obs}m_\lambda ^{mod}(E(BV),\tau )\}}{_\lambda 1/\sigma _\lambda }}.`$ (3)
with $`E(BV)`$ between 0.0 and 2.0 in steps of 0.02, where $`\sigma `$ is the uncertainty in the magnitude. This is a reasonable range, compared to the average colour excess for the bulge of 0.2 found by Lamers et al. (2002). With this initial guess we then make a three dimensional likelihood analysis, using a mass range from 1/1.5 to 1.5 times the initial mass estimate in steps of 0.004 dex. The extinction was varied between 0.0 and 2.0 in steps of 0.01 and the age was varied between 0.1 Myr and 1 Gyr for the Starburst99 models and between 10 Myr to 5 Gyr for the Frascati models.

For every fit we obtain a value for the reduced $`\chi ^2`$, i.e. $`\chi _\nu ^2`$ = $`\chi ^2`$/$`\nu `$, where $`\nu `$ is the number of free parameters i.e. the number of the observed data points minus the number of parameters in the theoretical model. For a good fit, $`\chi _\nu ^2`$ should be about unity. We checked that the fits were consistent with the faint magnitude limits of the filters in which the object was not detected. If not, the fit was rejected. The fit with the minimum value of $`\chi _\nu ^2`$ was adopted as the best fit and the corresponding values of $`E(BV)`$, age and $`M_i`$ were adopted. This method was applied for fits with the four sets of the Starburst99 models and with the Frascati models. Fig. 8 shows some examples of the results of the fitting process.

To estimate the uncertainty in the determined parameters we use confidence limits. If $`\chi _\nu ^2<\chi _\nu ^2(min)+1`$ then the resulting parameters, i.e. $`\mathrm{log}(t)`$, $`\mathrm{log}(M_i)`$ and $`E(BV)`$, are within the 68.3 % probability range. So the accepted ranges in age, mass and extinction are derived from the fits which have $`\chi _\nu ^2(min)<\chi _\nu ^2<\chi _\nu ^2(min)+1`$. With this method we derived the ages, initial masses333The initial cluster mass that we derived in this way is in fact “the initial mass of the current stellar population”, i.e. it is the “current mass corrected for stellar evolution effects”. This value can be different from the the initial mass of the original cluster if the cluster suffered evaporation or disruption. In that case the derived initial mass is the initial mass of the original cluster minus the fraction that has disappeared by evaporation or disruption. and extinction with their uncertainties of 602 clusters.

Fig. 9 shows a histogram of number of clusters as function of $`E(BV)`$. We will use this distribution for the two dimensional maximum likelihood fit of clusters detected in three bands only. The redding is small: 90 % of the clusters have a reddening lower than 0.40; 67 % of the clusters have a reddening lower than 0.18 and 23% have no detectable reddening.

### 5.3 Two dimensional maximum likelihood method

For the clusters observed in only three bands, it is not possible to make a three dimensional maximum likelihood fit. To reduce the parameter space we adopted the probability distribution of $`E(BV)`$ in the range of $`E(BV)`$ between 0.0 and 0.4, as shown in Fig. 9. For every value of $`E(BV)`$ between 0.0 and 0.40 (in steps of 0.02) and for every age of the model cluster, we first determine the mass of the cluster by means of Eq. (5.2). The results of the three dimensional fits have shown this is a good approximation. This results in a maximum likelihood age of the cluster for every value of $`E(BV)`$.

The $`\chi _\nu ^2`$ which comes out from the two dimensional maximum likelihood method is used to distinguish between the accepted and rejected fits. For the parameter $`\nu `$ we used a value of 3-1=2, because we only fit the age of the cluster. The mass of the cluster comes from the scaling of the magnitudes to the best-fit model. A fit is accepted if it agrees with the lower magnitude limits in the filters where it was not measured. To determine the age of the cluster, we average the ages, weighted by the probability that each particular value of $`E(BV)`$ occurs, derived by normalizing the distribution in Fig. 9 in the range of $`0.0<E(BV)<0.40`$. The error in the age is determined by calculating the value of the standard deviation $`\sigma `$ of the age, again weighted with the probability that the value of $`E(BV)`$ occurs. The value of $`E(BV)`$ is determined by using the one which belongs to the model with the age closest to the average age. The initial mass of the cluster is then calculated by Eq. (5.2) for the adopted value of $`E(BV)`$. With this method we derived the age, initial mass and the extinction of 275 clusters, using both the Starburst99 models and the Frascati models.

### 5.4 Starburst99 or Frascati models?

We have compared the energy distributions of the objects with five sets of models: four sets of Starburst99 models, with $`Z=Z_{}`$ or $`2Z_{}`$ and $`M_{\mathrm{up}}=30`$ or $`100`$ $`M_{}`$, and the Frascati models with $`Z=Z_{}`$ and $`M_{\mathrm{up}}=25`$ $`M_{}`$. Based on the lack of \[OIII\] emission we adopt the models with $`M_{\mathrm{up}}=25`$ or 30 $`M_{}`$. The models with $`M_{\mathrm{up}}=100`$ $`M_{}`$ are only used later to check the influence of this adopted upper limit on our results.

We have compared the ages derived from fitting the observed SEDs with the Starburst99 models and the Frascati models, for objects measured in at least four filters and with fits of $`\chi _\nu ^2<3.0`$. We found that in the age range of 10 to 800 Myr there is a reasonable correlation between the results of the SB99 models and the Frascati models, with the Frascati models giving ages systematically about 0.4 dex smaller than the SB99 models. In this age range we adopt the SB99 models, because they are based on more reliable evolutionary tracks, better stellar atmosphere models and because the nebular continuum is included. We find that in this mass range the fits with the SB99 models have smaller $`\chi _\nu ^2`$ than the fits with the Frascati-models. For objects with ages above about 700 Myrs ($`\mathrm{log}(t)>8.85`$) the Frascati-models give the most reliable fits with smaller $`\chi _\nu ^2`$ than the SB99 models. This is probably because the lower mass limit of the SB99 models is 1 $`M_{}`$, and the Frascati models go down to 0.6 $`M_{}`$. Moreover the SB99 models do not go beyond 1 Gyr. Based on this comparison we adopt the SB99 fits for ages less than 700 Myrs and the results of the Frascati models for older ages.

### 5.5 Uncertainties in the derived parameters

The determination of the cluster parameters, age, mass and extinction, by means of the two or thee dimensional maximum likelihood fitting method, i.e. the 3(2)DEF-method, results not only in the values of $`\mathrm{log}(M_{\mathrm{cl}})`$, $`\mathrm{log}(t)`$ and $`E(BV)`$ but also in their maximum and minimum acceptable values. The best-fit values are not necessarily in the middle of the minimum and the maximum values. If we define the uncertainties in the parameters, $`\mathrm{\Delta }E(BV)`$, $`\mathrm{\Delta }\mathrm{log}(M_{\mathrm{cl}})`$ and $`\mathrm{\Delta }\mathrm{log}(t)`$ as half the difference between the maximum and mimimum values, we find that 30 percent of the clusters have $`\mathrm{\Delta }E(BV)<0.08`$, 50 percent have $`\mathrm{\Delta }E(BV)<0.11`$ and 70 percent have $`\mathrm{\Delta }E(BV)<0.15`$. For the uncertainties in the mass determination we find $`\mathrm{\Delta }\mathrm{log}(M_{\mathrm{cl}})<0.20`$, $`<0.33`$ and $`<0.55`$ for 30, 50 and 70 percent of the clusters respectively. For the uncertainties in the age determination we find $`\mathrm{\Delta }\mathrm{log}(t)<0.23`$, $`<0.39`$ and $`<0.69`$ for 30, 50 and 70 percent of the clusters respectively. Taking the values for 50 percent of the clusters as representative, we conclude that the uncertainties are $`\mathrm{\Delta }E(BV)0.11`$, $`\mathrm{\Delta }\mathrm{log}(M_{\mathrm{cl}})0.33`$ and $`\mathrm{\Delta }\mathrm{log}(t)0.39`$.

### 5.6 Contamination of the cluster sample by stars?

To estimate the number of stars that may contaminate our cluster sample we use the stellar population in the solar neighbourhood. ¿From the tabulated stellar densities as a function of spectral type (Allen 1976) we derived the number of stars brighter than a certain value of $`M_V`$ per pc3. Using the mass density of 0.13 $`M_{}`$pc-3, we derived the number of stars per unit stellar mass. The results are listed in Table 2. The total stellar mass of the observed region of M51 (see Fig. 2) is estimated to be about $`1/20`$ of the total mass of $`\mathrm{5\; 10}^{10}`$ $`M_{}`$ in the disk of that galaxy (Athanassoula et al. 1987), i.e. about $`\mathrm{2.5\; 10}^9`$ $`M_{}`$. This implies that we can expect the following numbers of stars brighter than $`M_V<6.5`$ per spectral type in the observed region: 4 (O), 13 (B), 6 (A), 6 (F), 6 (G), $`<`$1 (K) and $`<`$1 (M). So we can expect of the order of 40 bright stars with $`M_V<6.5`$ in the observed region of M51. However, these are all massive young supergiants of which the vast majority will be in clusters! So the number of bright stars outside clusters, that may contaminate our sample of clusters will be considerable smaller, and we expect it to be smaller than about 20 out of the total sample of 877. Moreover, we have shown above, from the lack of $`O[III]`$ emission, that the clusters in M51 have a shortage of massive stars with $`M>30`$ $`M_{}`$. We can expect this effect also to occur for the massive field stars. This would reduce the the number of expected contaminating stars even further. Tests have shown that a considerable fraction of possibly remaining stellar sources will be eliminated by the requirement that their energy distribution should fit that of cluster models within a given accuracy. Based on all these considerations, we conclude that contamination of our cluster sample with very massive stars of $`M_V<6.5`$ outside clusters is expected to be negligible.

## 6 The mass versus age distribution

We discuss the results from the fitting of the observed energy distributions to those of cluster models with solar metallicity and with and upper mass limit of 30 $`M_{}`$. We have checked that the fits with twice solar metallicity and with an upper limit of 100 $`M_{}`$ give about the same results.

Fig. 10 shows the mass versus age distribution of sources with $`M_v<7.5`$ (to eliminate possible stellar sources) and with energy distributions that could be fitted to that of a cluster models with an accuracy of $`\chi _\nu ^23.0`$. (The distributions for clusters with $`\chi _\nu ^21`$ or 10, not shown here, show the same distribution, but with 294 and 508 clusters respectively.) We see that the lower limit of the mass increases with increasing age, from about 1000 $`M_{}`$ at $`\tau 5`$ Myrs to about $`\mathrm{5\; 10}^4`$ $`M_{}`$ at 1 Gyr. This is due to the expected effect of fading of the clusters as they age (see Sect. 4.1). The full line in the figure is the fading line in the $`R`$ magnitude for clusters which have a limiting magnitude of $`R_{\mathrm{lim}}=22.0`$ at the distance of M51 for a reddening of $`E(BV)`$=0. The line in Fig. 10 thus gives the mass of the clusters, that reaches this magnitude detection limit, as a function of age.

These initial mass versus age distribution of the clusters in Fig. 10 show the following characteristics:
– (i) The lower mass limit increases with age due to the fading of the clusters. The observed lower limit agrees with the predicted ones for the $`R`$ band. There are even hints of the presence of the predicted dips in the lower limit near $`\mathrm{log}(t)`$=6.8 and 7.2. This strengthens our confidence in the adopted models.
– (ii) There are clear concentrations in the distribution at $`\mathrm{log}(t)=6.70`$ and 7.45 and possibly also around $`\mathrm{log}(t)7.2`$. These are due to the properties of the cluster models and the adopted method. In Fig. 7 we have shown that the colours of the models do not change monotonically with time, but that there are phases when the colours change rapidly with time. These phases occur in the range of $`6.5<\mathrm{log}(t)<7.3`$. Large changes in the colours of the models occur just before and after the peaks where the slopes of the curves in Fig. 7 are large. It is more difficult to fit the observed energy distributions with high accuracy to models in the age range when the spectral changes are large, than in the age range when the changes are small. So there is a tendency of the model fits to concentrate in agebins just outside the age-regions of large spectral changes. This explains the concentrations and the voids in the derived age distributions of the clusters between about $`6.5<\mathrm{log}(t)<7.5`$.
– (iii) The density of the points drops at $`\mathrm{log}(t)>7.5`$. Since the age scale is logarithmic, we might have expected an increasing density of points towards the higher ages, which is not observed. This is due to the disruption or dispersion of the clusters (see first paragraphs of Sect. 7).

We conclude that the mass versus age distribution agrees with the expected evolutionary fading of the clusters and that the decrease in numbers of clusters with age shows the affect of disruption/dispersion of clusters with time. The concentrations of the clusters at ages around $`\mathrm{log}(t)=6.8`$, 7.2 and 7.45 are due to statistical effects and do not represent periods of enhanced cluster formation.

## 7 The cluster Initial Mass Function

The data in Fig. 10 and the detailed study of this distribution by Boutloukos & Lamers (2002) show that the clusters in the inner spiral arms of M51 disrupt on a time scale of about tens of Myrs. In fact, these authors derived the dependence of the disruption time on the initial mass of the clusters in the inner spiral arms of M51 as

$$\mathrm{log}t_{\mathrm{disr}}=\mathrm{log}t_4+\gamma \times \mathrm{log}(M_{\mathrm{cl}}/10^4M_{})$$
(4)
with $`\mathrm{log}t_4=7.64\pm 0.22`$ and $`\gamma =0.62\pm 0.06`$ for the mass range of $`3\mathrm{log}(M_{\mathrm{cl}}/M_{})5.2`$, where $`M_{\mathrm{cl}}`$ is the initial mass of the cluster. We see that clusters with an initial mass larger than $`10^4`$ $`M_{}`$ survive $`4\times 10^7`$ years. Clusters with an initial mass of only $`10^3`$ $`M_{}`$ disrupt on a time scale of $`1\times 10^7`$ yrs. This implies that the cluster initial mass function cannot be derived from the total sample of clusters, because the disruption will produce a strong bias towards the more massive clusters. However for the youngest clusters with ages less than about 10 Myr disruption is not yet an important effect and, therefore, these clusters can be used to derive the initial cluster mass function.

Fig. 11 shows the resulting mass distribution of clusters with an age less or equal to 10 Myr for two samples of clusters. The first sample contains all 354 clusters younger than 10 Myr. The second sample contains 168 clusters in the same age range but with a mass determination of $`\mathrm{\Delta }\mathrm{log}(M_{\mathrm{cl}})0.25`$. Both samples show the same characteristics: a steep increase in number between $`2.5<\mathrm{log}M_{\mathrm{cl}}<3.0`$ and a slow decrease to higher masses. The steep increase is due to the detection limit or the disruption of the low mass clusters. The slow decrease reflects the cluster initial mass function (CIMF). The decrease indicates that the clusters are formed with an CIMF that has a negative slope of $`d\mathrm{log}(N)/d\mathrm{log}(M)`$, as expected.

If the CIMF can be written as a power law of the type

$$N(M)dMM^\alpha dM\mathrm{for}M_{\mathrm{min}}<M<M_{\mathrm{max}}$$
(5)
then the normalized cumulative distribution will be

$`\mathrm{\Sigma }(M)`$ $`=`$ $`AB\times M^{\alpha +1}\mathrm{with}`$
$`A`$ $`=`$ $`B\times M_{\mathrm{min}}^{\alpha +1}`$
$`B`$ $`=`$ $`(M_{\mathrm{min}}^{\alpha +1}M_{\mathrm{max}}^{\alpha +1})^1`$ (6)
Fig. 12 shows the observed and predicted normalized cumulative distribution of the 149 clusters with an age less than 10 Myr and with a mass $`M_{\mathrm{cl}}>\mathrm{2.5\; 10}^3M_{}`$ that are used for the determination of the CIMF. We eliminated the clusters with $`M_{\mathrm{cl}}<\mathrm{2.5\; 10}^3M_{}`$ from the sample because Fig. 11 shows that the sample may not be complete for smaller masses. The cumulative distributions of the 84 clusters of $`M>\mathrm{2.5\; 10}^3M_{}`$, younger than 10 Myr, which are fitted to the models with an accuracy of $`\chi _\nu ^23.0`$, or of the 66 clusters with $`\mathrm{\Delta }\mathrm{log}(M_{\mathrm{cl}})0.25`$, not shown here, (where $`\mathrm{\Delta }\mathrm{log}(M_{\mathrm{cl}})`$ is the uncertainty in the mass determination), have the same shape as the distribution in Fig. 12. For these three samples we have determined the values of $`\alpha `$ and $`M_{\mathrm{min}}`$ by means of a linear regression under the reasonable assumption that $`M_{\mathrm{min}}<<M_{\mathrm{max}}`$, which implies that $`A1`$. The resulting values of $`\alpha `$ and $`M_{\mathrm{min}}`$ are listed in Table 3. The table shows that $`\alpha 2.1`$ for all three samples.

The predicted cumulative distribution with the derived values of $`\alpha =2.12`$ and $`M_{\mathrm{min}}`$ is shown in Fig. 12 (full line). The figure shows a slight underabundance of clusters in the range of $`4.0<\mathrm{log}(M_{\mathrm{cl}})<4.3`$ and a slight overabundance in the range of $`\mathrm{log}(M_{\mathrm{cl}})>4.3`$. In fact, for the mass range of $`3.5<\mathrm{log}(M_{\mathrm{cl}})<4.3`$ a fit with $`\alpha =2.00`$, shown by dashed lines, fits the distribution excellently. We conclude that the CIMF of clusters younger than 10 Myr has a slope of $`\alpha 2.1\pm 0.3`$ in the mass range of $`3.0<\mathrm{log}(M_{\mathrm{cl}})<5.0`$ and a slope of $`2.00\pm 0.05`$ in the range of $`3.0<\mathrm{log}(M_{\mathrm{cl}})<4.3M_{}`$.

The derived exponent of the cluster IMF is very similar to the value of $`\alpha =2.0`$ of young clusters in the Antennae galaxies, as found by Zhang & Fall (1999). It shows that the IMF of clusters formed in the process of galaxy-galaxy interaction is very similar to the one of clusters formed in the spiral arms of a galaxy, long after the interaction. This mass distribution is also similar to that of giant molecular clouds (e.g. McKee 1999; Myers 1999). This may support the suggestion that the mass distribution of the clusters is determined by the mass distribution of the clouds from which they originate.

## 8 The cluster formation history

One of the goals of this paper is to study the influence of the interaction between M51 and its companion on the cluster formation in the region of the inner spiral arms. To this purpose we compare the observed age distribution of the clusters with predictions for a constant cluster formation rate. This comparison is hampered by two effects: (a) the disruption of clusters and (b) the fading of clusters below the detection limit. The disruption time of the clusters in the inner spiral arms of M51 is given by Eq. (4). We see that only clusters with an initial mass larger than about $`10^4`$ $`M_{}`$ will survive more than 40 Myrs.

Figure 13 shows the histogram of the formation rates of the detected clusters, in number per Myr, for 140 clusters with $`M_{\mathrm{cl}}>10^4`$ $`M_{}`$ and with an energy distribution that is fitted with an accuracy of $`\chi _\nu ^23.0`$. The sample of 111 clusters with $`M_{\mathrm{cl}}>10^4`$ $`M_{}`$ and with an age uncertainty of $`\mathrm{\Delta }\mathrm{log}(t)<0.25`$, and the full sample of 285 clusters with $`M_{\mathrm{cl}}>10^4`$ $`M_{}`$, not shown here, have a distribution very similar to the one shown in Fig. 12.

All three samples show about the same characteristics.
(a) There is a general trend of a decrease in the cluster formation rate towards increasing age.
(b) There is a steep drop around $`\mathrm{log}(t)7.5`$. This drop is the result of the concentration of clusters at ages $`\mathrm{log}(t)=7.45`$ that was apparent in the mass versus age diagrams of Fig. 10, and was explained in Sect. 6.
(c) There is no clear evidence for a peak in the cluster formation rate near $`\mathrm{log}(t)8.6`$, which is the time of the interaction of M51 with its companion, and which is the age of the huge starburst in the nucleus of M51. There is a hint of a small peak around $`\mathrm{log}(t)8.1`$. However this peak is only $`2\sigma `$ high. Its reality has to be verified with a larger sample of clusters.
(d) The peak in the last bin of 5 Gyr contains all the clusters ages older than 3 Gyr, because the cluster models that we used for fitting the energy distributions do not go beyond 5 Gyr.
We have checked that these characteristics are not the result of the binning process: the same features appear for different choices of the binning parameters.

The general decrease in the formation rate of the observed clusters is partly due to the evolutionary fading of the clusters and partly due to the disruption of clusters. In Sect. 6 and in Fig. 10 we have shown that clusters with an initial mass of $`10^4`$ $`M_{}`$ fade below the detection limit when they are older than $`\mathrm{log}(t)=8.2`$. This can explain the decrease in the formation rate of clusters older than about 150 Myrs. However, clusters with $`M_{\mathrm{cl}}>10^4M_{}`$ and younger than about 100 Myr should still all be detectable. The fact that we see a decrease in the apparent cluster formation rate must thus be due to the disruption of clusters (unless for some unknown reason the cluster formation rate has been steadily increasing from 100 Myrs up to now, which we consider unlikely).

It is interesting that the CFR in the inner spiral arms of M51, at a distance of about 1 to 3 kpc from the nucleus, does not show any evidence for a peak at the age of the strongest interaction of M51 with its companion galaxy. The closest approach occurred about 250 - 400 Myrs ago according to the dynamical models of Barnes (1998) or 400-500 and 50-100 Myrs ago according to Salo & Laurikainen (2000). The only significant peak occurred at $`\mathrm{log}(t)7.4`$. However we attributed this peak to the large changes in the energy distributions of the models in the age range of $`6.5<\mathrm{log}(t)<7.5`$, and not to a real increase in the CFR.

## 9 Summary and conclusions

We have studied images of the inner spiral arms of the interacting galaxy M51 obtained with the $`HSTWFPC2`$ camera in broad-band $`UBVRI`$ and in the narrow band H$`\alpha `$ and \[OIII\] filters. The study can be summarized as follows:

1. We found a total of 877 point-like objects, which are probably clusters. Many of the clusters are strong H$`\alpha `$ emitters, but none of the clusters, not even the youngest ones, have an excess of radiation in the \[OIII\] line at 5007 Å ($`F502N`$-filter). This suggests that the upper mass limit of the stars in the clusters is about 25 to 30 $`M_{}`$.
2. We have compared their energy distributions with those of Starburst99 cluster models (Leitherer et al. 1999) for instantaneous star formation with a stellar IMF of exponent 2.35, solar metallicity, a lower and upper stellar mass limit of 1 $`M_{}`$ and 30 $`M_{}`$ respectively. The energy distributions were also compared with those of the Frascati models (Romaniello 1998). For clusters younger than 700 Myr the results from the fitting with the Starburst99 models were adopted because these models are more accurate for young clusters and the fits of the energy distributions are better than those of the Frascati models. For older clusters the results from the fits with the Frascati models were adopted.
3. For clusters that were observed in four or five bands a three dimensional maximum likelihood method was used to derive the properties of the clusters from the comparison between the observed and predicted energy distributions. The free parameters are the age $`t`$ and $`E(BV)`$, which together determined the shape of the energy distribution, and the initial cluster mass $`M_{\mathrm{cl}}`$ which determines the absolute magnitude. For clusters that were not observed in all bands, the empirically derived lower magnitudes limits were taken into account.
4. For clusters that were observed in only three bands the age and mass were derived in a two-dimensional maximum likelihood fitting of the energy distributions, with $`t`$ and $`M_{\mathrm{cl}}`$ as free parameters. The observed probability distribution of $`E(BV)`$ was used as a weighting factor in the fitting procedure.
5. The histogram of $`E(BV)`$ is strongly peaked at very small $`E(BV)0`$. All cluster have a reddening smaller than $`E(BV)<1.0`$ and 67 % of the clusters have $`E(BV)<0.18`$.
6. We have analysed the observed clusters also with cluster models of higher metallicity, $`Z=2\times Z_{}`$. These higher metallicity models fit the observations considerably worse than the solar metallicity models. For instance, for solar metallicity models the energy distribution of 294 clusters can be fitted with an accuracy of $`\chi _\nu ^21.0`$ and 392 with $`\chi _\nu ^23.0`$. For models with twice the solar metallicity these numbers are respectively 138 and 217. So the energy distributions of the clusters support the adopted solar metallicity.
7. The clusters have masses in the range of $`2.5<\mathrm{log}(M_{\mathrm{cl}})<5.7`$ and ages of $`\mathrm{log}(t)>5.0`$. These masses are the initial masses of the clusters, i.e. the current mass corrected for stellar evolution effects, but not corrected for evaporation or disruption. All derived masses have to be multiplied by a factor 1.3 if the lower mass of the stars is 0.6 $`M_{}`$, instead of the adopted 1 $`M_{}`$, and by a factor 2.1 if the lower mass is 0.2 $`M_{}`$, as found for the Orion Nebula cluster.
8. The distribution of the clusters in a mass-versus-age diagram shows the predicted lower limit due to the evolutionary fading of the clusters, including the dips at $`\mathrm{log}(t)6.8`$ and 7.1. Three apparent concentrations at $`\mathrm{log}(t)=6.7`$, 7.2 and 7.45 are not real but due to the properties of the cluster models used.
9. About 60% of the clusters are younger than 40 Myr. The number of older clusters is much less than expected for a constant cluster formation rate. This is partly due to the evolutionary fading of low mass clusters below the detection limit, and partly due to the disruption of the clusters.
10. The cluster initial mass function (CIMF) was derived from the cumulative mass distribution of clusters younger than 10 Myr, for which disruption has not occured. The CIMF has a slope of $`\alpha =2.1\pm 0.3`$ in the range of $`3.0<\mathrm{log}(M_{\mathrm{cl}})<5.0`$ and $`\alpha =2.00\pm 0.05`$ in the range of $`3.0<\mathrm{log}(M_{\mathrm{cl}})<4.5`$ $`M_{}`$, for $`N(M_{\mathrm{cl}})M_{\mathrm{cl}}^\alpha `$. This slope is the same to that found in the interacting Antennae galaxies (Zhang & Fall 1999). Zhang and Fall deived a power law slope of the CIMF of $`\alpha =1.95\pm 0.03`$ and $`2.00\pm 0.04`$ for two cluster samples of the Antennae galaxies. The good agreement between these slopes and the one found by us suggests that $`\alpha `$ is about the same for cluster formation triggered by strong galaxy-galaxy interactions, such as presently going on in the Antennae, as for cluster formation that is not dominated by the interactions.
11. The age distribution of clusters with $`M_{\mathrm{cl}}>10^4`$ $`M_{}`$, is used to derive the history of the cluster formation rate (CFR). There is a general trend of a decrease of the formation rate of the observed clusters with age. It is unlikely that the real CFR has been increasing continuously from about 1 Gyr to the present time. The decrease of the CFR with age of clusters younger than about 100 Myr cannot be due to evolutionary fading, but it is due to the disruption of clusters. For clusters older than 200 Myr the decrease of the derived CFR could, at least partly, be due to evolutionary fading.
12. There is no evidence for a peak in the CFR at about 400 Myr, which is the time of the interaction of M51 with its companion and the age of the huge starburst in the nucleus.

In a forthcoming paper we describe the cluster formation as a function of location in a large part of M51, using the same methods as used here (Bastian et al. 2002). The disruption of clusters in M51, derived from the results of the study presented here, are described by Boutloukos & Lamers (2002).

## 10 Acknowledgement

H.J.G.L.M.L. and N.B. are grateful to the Space Telescope Scence Institute for hospitality and financial support during several stays. We thank Claus Leitherer for help and advice in the calculation of the cluster models. Support for the SINS program GO-9114 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. N.B. ackowledges a grant from the Netherlands Organization for Scientific Research. We thank the unknown referee for constructive comments that resulted in an improvement of this paper.