text:
# Ab Initio Approach to 𝑠-Shell Hypernuclei ³{_Λ}H, ⁴{_Λ}H, ⁴{_Λ}He and ⁵{_Λ}He with a Λ𝑁-Σ𝑁 Interaction

## Abstract

Variational calculations for $`s`$-shell hypernuclei are performed by explicitly including $`\mathrm{\Sigma }`$ degrees of freedom. Four sets of $`YN`$ interactions (SC97d(S), SC97e(S), SC97f(S) and SC89(S)) are used. The bound-state solution of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He is obtained and a large energy expectation value of the tensor $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ transition part is found. The internal energy of the 4He subsystem is strongly affected by the presence of a $`\mathrm{\Lambda }`$ particle with the strong tensor $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ transition potential.

Few-body calculations for $`s`$-shell hypernuclei with mass number $`A=35`$ are important not only to explore exotic nuclear structure, including the strangeness degrees of freedom, but also to clarify the characteristic features of the hyperon-nucleon ($`YN`$) interaction. Although several interaction models have been proposedNSC89 ; NSC97 ; FSS96 , the detailed properties (e.g. $`{}_{}{}^{1}S_{0}^{}`$ or $`{}_{}{}^{3}S_{1}^{}^3D_1`$ phase shift, strength of $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ coupling term) of the $`YN`$ interaction are different among the models. The observed separation energies ($`B_\mathrm{\Lambda }`$’s) of light $`\mathrm{\Lambda }`$ hypernuclei are expected to provide important information on the $`YN`$ interaction, because the relative strength of the spin-dependent term or of the $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ coupling term is affected from system to system.

Recently, few-body studies for $`A=3,4`$ hypernuclei have been conducted using modern $`YN`$ interactionsMiyagawa ; Hiyama ; Nogga . According to these developments, the Nijmegen soft core (NSC) model 97f (or 97e) seems to be compatible with the experimental $`B_\mathrm{\Lambda }`$’s, though the calculated $`B_\mathrm{\Lambda }`$ for $`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$H or $`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$He is actually slightly smaller than the experimental value. These few-body calculations, however, have not yet reached a stage to calculate $`B_\mathrm{\Lambda }`$($`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He).

If one constructs a phenomenological central $`\mathrm{\Lambda }N`$ potential, which is consistent with the experimental $`B_\mathrm{\Lambda }`$($`{}_{}{}^{3}{}_{\mathrm{\Lambda }}{}^{}`$H), $`B_\mathrm{\Lambda }`$($`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$H), $`B_\mathrm{\Lambda }`$($`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$He), $`B_\mathrm{\Lambda }`$($`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$H) and $`B_\mathrm{\Lambda }`$($`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$He) values as well as the $`\mathrm{\Lambda }p`$ total cross section, that kind of potential would overestimate the $`B_\mathrm{\Lambda }`$($`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He) valueDHT ; Nemura . This is known as an anomalously small binding of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He. Though a suppression of the tensor forces Bando ; Shinmura  or of the $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ coupling Bodmer ; Gibson  was discussed to be a possible mechanism to resolve the anomaly, the problem still remains an enigmaSURVEY  due to the difficulty of performing a complete five-body treatment. Only one attempt was made, using a variational Monte Carlo calculationCarlson  with the NSC89 $`YN`$ interaction. Though NSC89 well reproduces both the experimental $`B_\mathrm{\Lambda }`$($`{}_{}{}^{3}{}_{\mathrm{\Lambda }}{}^{}`$H)Miyagawa  and $`B_\mathrm{\Lambda }`$($`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$H)Carlson ; Nogga  values as well as the experimental $`\mathrm{\Lambda }p`$ total cross section, a bound-state solution of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He was not found. In view of the aim to pin down a reliable $`YN`$ interaction, a systematic study for all $`s`$-shell hypernuclei is desirable.

The $`NN`$ tensor interaction due to a one-pion-exchange mechanism is the most important ingredient for the binding mechanisms of light nuclei. More than a third, or about one half, of the interaction energy comes from the tensor force for the 4HeATMS ; SVM ; Pieper . Since the pion-(or kaon-) exchange also induces the $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ transition for the $`YN`$ sector, both the $`NN`$ and $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ tensor interactions may also play important roles for light hypernuclei. If this is the case, the structure of the core nucleus (e.g. 4He) in the hypernucleus ($`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He) would be strongly influenced by the presence of a $`\mathrm{\Lambda }`$ particle.

The purpose of this letter is twofold: First is to perform an ab initio calculation for $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He as well as $`A=3,4`$ hypernuclei explicitly including $`\mathrm{\Sigma }`$ degrees of freedom. Second is to discuss the structural aspects of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He with an appropriate $`YN`$ interaction which is consistent with all of the $`s`$-shell hypernuclear data.

The Hamiltonian ($`H`$) of a system comprising nucleons and a hyperon ($`\mathrm{\Lambda }`$ or $`\mathrm{\Sigma }`$) is given by $`2\times 2`$ components as

$$H=\left(\begin{array}{cc}H_\mathrm{\Lambda }& V_{\mathrm{\Sigma }\mathrm{\Lambda }}\\ V_{\mathrm{\Lambda }\mathrm{\Sigma }}& H_\mathrm{\Sigma }\end{array}\right),$$
(1)
where $`H_\mathrm{\Lambda }(H_\mathrm{\Sigma })`$ operates on the $`\mathrm{\Lambda }`$-($`\mathrm{\Sigma }`$-)component and

$$V_{\mathrm{\Lambda }\mathrm{\Sigma }}=\underset{i=1}{\overset{A1}{}}v_{iY}^{(N\mathrm{\Lambda }N\mathrm{\Sigma })}.$$
(2)
We employ the G3RS potentialTamagaki  for the $`NN`$ interaction and the SC97d(S), SC97e(S), SC97f(S) or SC89(S) potentialAkaishi  for the $`YN`$ interaction, where all interactions have tensor and spin-orbit components in addition to the central one. We omit small nonstatic correction terms ($`(𝑳𝑺)^2`$ and $`𝑳^2`$ terms) in the G3RS $`NN`$ interaction and odd partial-wave components in each interaction in order to focus on the main part of the interaction in the even parity state. The calculated binding energies for light nuclei (2H, 3H, 3He and 4He) are 2.28, 7.63, 6.98 and 24.57 MeV, respectively. The $`YN`$ interactions have Gaussian form factors whose parameters are set to reproduce the low-energy $`S`$ matrix of the corresponding original Nijmegen $`YN`$ interactionsShinPrivate . These Gaussian form factors help to save significant computer time.

The binding energies of various systems are calculated in a complete $`A`$-body treatment. The variational trial function must be flexible enough to incorporate both the explicit $`\mathrm{\Sigma }`$ degrees of freedom and higher orbital angular momenta. The trial function is given by a combination of basis functions:

$$\begin{array}{c}\mathrm{\Psi }_{JMTM_T}=_{k=1}^Nc_k\phi _k,\text{with}\hfill \\ \phi _k=𝒜\left\{G(𝒙;A_k)[\theta _{L_k}(𝒙;u_k,K_k)\chi _{S_k}]_{JM}\eta _{kTM_T}\right\}.\hfill \end{array}$$
(3)
Here, $`𝒜`$ is an antisymmetrizer acting on nucleons and $`\chi _{S_k}`$ $`\left(\eta _{kTM_T}\right)`$ is the spin (isospin) function. $`\eta _{kTM_T}`$ has two components: upper (lower) component refers to the $`\mathrm{\Lambda }`$-($`\mathrm{\Sigma }`$-)component. The abbreviation $`𝒙=(𝒙_1,\mathrm{},𝒙_{A1})`$ is a set of relative coordinates. A set of linear variational parameters $`(c_1,\mathrm{},c_N)`$ is determined by the Ritz variational principle.

A spatial part of the basis function is constructed by the correlated Gaussian(CG) multiplied by the orbital angular momentum part $`\theta _L(𝒙)`$, expressed by the global vector representation(GVR)GlobalV . CG is defined by

$$\begin{array}{cc}\hfill G(𝒙;A_k)& =\mathrm{exp}\left\{\frac{1}{2}_{i<j}^A\alpha _{k}^{}{}_{ij}{}^{}(𝒓_i𝒓_j)^2\right\}\hfill \\ & =\mathrm{exp}\left\{\frac{1}{2}_{i,j=1}^{A1}(A_k)_{ij}𝒙_i𝒙_j\right\}.\hfill \end{array}$$
(4)
The $`(A1)\times (A1)`$ symmetric matrix ($`A_k`$) is uniquely determined in terms of the interparticle correlation parameter ($`\alpha _{kij}`$). The GVR of $`\theta _{L_k}(𝒙;u_k,K_k)`$ takes the form

$$\begin{array}{c}\theta _{L_k}(𝒙;u_k,K_k)=v_k^{2K_k+L_k}Y_{L_k}(\widehat{𝒗}_k),\text{with}\hfill \\ 𝒗_k=_{i=1}^{A1}(u_k)_i𝒙_i.\hfill \end{array}$$
(5)
The $`A_k`$ and $`u_k`$ are sets of nonlinear parameters which characterize the spatial part of the basis function. Allowing the factor $`v_k^{2K_k}`$ $`(K_k0)`$ is useful to improve the short-range behavior of the trial function. The value of $`K_k`$ is assumed to take 0 or 1. The variational parameters are optimized by a stochastic procedure. The above form of the trial function gives accurate solutions. The reader is referred to Ref.SVM ; GlobalV  for details and recent applications. For the spin and isospin parts, all possible configurations are taken into account.

Table 1 lists the results of the $`\mathrm{\Lambda }`$ separation energies. The scattering lengths of the $`{}_{}{}^{1}S_{0}^{}(a_s)`$ and $`{}_{}{}^{3}S_{1}^{}(a_t)`$ states for each $`YN`$ interaction are also listed in Table 1, where the interactions are given in increasing order of $`|a_s|`$ (and in decreasing order of $`|a_t|`$). The SC89(S) interaction produces no or very weakly bound state for $`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$H, $`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$He or $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He. For the SC97d-f(S), the $`B_\mathrm{\Lambda }`$($`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He) value is about $`23`$ MeV. This is a first ab initio calculation to produce the bound state of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He with explicit $`\mathrm{\Sigma }`$ degrees of freedom.

The order of the spin doublet structure of the $`A=4`$ system is correctly reproduced for all $`YN`$ interactions; the ground (excited) state has spin-parity, $`J^\pi =0^+(1^+)`$ for both isodoublet hypernuclei $`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$H and $`{}_{}{}^{4}{}_{\mathrm{\Lambda }}{}^{}`$He. Although the strengths of the $`{}_{}{}^{1}S_{0}^{}`$ and $`{}_{}{}^{3}S_{1}^{}`$ interactions of the SC97d(S) are almost the same as each other, the energy-level of the $`0^+`$ state is clearly lower than that of the $`1^+`$ state. All of the $`A=3`$ bound states given in Table 1 have $`J^\pi =\frac{1}{2}^+`$, in agreement with experiment. No other bound state has been obtained for all of the $`YN`$ interactions. For the SC97e(S), the differences between the calculated and experimental $`B_\mathrm{\Lambda }`$ values are the smallest among the $`YN`$ interactions employed in the present study.

Table 2 lists the probability, $`P_\mathrm{\Sigma }`$ (in percentage), of finding a $`\mathrm{\Sigma }`$ particle in the system. The sizable amount of $`P_\mathrm{\Sigma }`$($`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He)’s is obtained. This implies that the $`\mathrm{\Lambda }\mathrm{\Sigma }`$ coupling plays an important role, even for the $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He, despite a large excitation energy of the core nucleus, 4He (with the isospin 1), in the $`\mathrm{\Sigma }`$-component. For the $`A=4`$ system, the $`P_\mathrm{\Sigma }`$’s of the $`0^+`$ state are about $`12\%`$, except for the SC89(S), while the $`P_\mathrm{\Sigma }`$’s of the $`1^+`$ state are nearly equal to or smaller than that of the $`0^+`$ state.

Figure 1 displays the density distributions for $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He using SC97e(S), and of $`N,\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$ from the center-of-mass (CM) of 4He. Figure 1 also shows the $`\mathrm{\Lambda }`$-distribution obtained from the Isle $`\mathrm{\Lambda }\alpha `$ potentialFuse . The experimental pionic decay width of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He suggests that the $`\mathrm{\Lambda }`$-distribution should spread over a rather outer region compared to the distribution of the $`\alpha `$, as was discussed in Ref. Fuse . The present curve of the $`\mathrm{\Lambda }`$-distribution is similar to that obtained by the Isle potential. The $`\mathrm{\Sigma }`$-distribution has a shape similar to the $`N`$-distribution. The root-mean-square (rms) radii of $`N`$, $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$ from the CM of the 4He are $`1.5`$, $`2.9`$ and $`1.6`$ fm, respectively.

Table 3 lists the energy expectation values of the kinetic and potential energy terms for $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He. The contributions from the spin-orbit and the Coulomb potentials are not shown in the table, though the calculations include them. Here, $`T_c`$ is the kinetic energy of the core nucleus ($`c`$) subtracted by the CM energy of $`c`$:

$$T_c=\underset{i=1}{\overset{A1}{}}\frac{𝒑_i^2}{2m_N}\frac{\left(\underset{i=1}{\overset{A1}{}}𝒑_i\right)^2}{2(A1)m_N}.$$
(6)
The kinetic energy of the relative motion between the $`Y`$ and the CM of $`c`$ is given by

$$T_{Yc}=\frac{𝝅_{Yc}^2}{2\mu _Y}+(m_Ym_\mathrm{\Lambda })c^2,$$
(7)
where $`\mu _Y=\frac{(A1)m_Nm_Y}{(A1)m_N+m_Y}`$ is the reduced mass for the $`Y+c`$ system and $`𝝅_{Yc}`$ is the canonical momentum of the relative coordinate between $`Y`$ and $`c`$ $`(Y=\mathrm{\Lambda },\mathrm{\Sigma })`$. $`T_{Yc}`$ also counts the difference in the rest-mass energy between $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$. Each potential part $`V`$ takes account of a summation over appropriate particle pairs (see Eq. (2) for example). The energy expectation values of the first three columns in Table 3 are written as

$$\begin{array}{ccc}\hfill 𝒪& =& \mathrm{\Psi }_\mathrm{\Lambda }|𝒪|\mathrm{\Psi }_\mathrm{\Lambda }+\mathrm{\Psi }_\mathrm{\Sigma }|𝒪|\mathrm{\Psi }_\mathrm{\Sigma },\hfill \end{array}$$
(8)
where the upper (lower) component of the $`\mathrm{\Psi }_{JMTM_T}`$ is denoted by $`\mathrm{\Psi }_\mathrm{\Lambda }`$ ($`\mathrm{\Psi }_\mathrm{\Sigma }`$). The first (second) term of each element ($`T_c`$, $`T_{Yc}`$ or $`V_{NN}`$) in Table 3 represents the first (second) term of Eq. (8). The energy of the 4He subsystem changes a lot from that of the isolated one,

$`\mathrm{\Delta }E_c`$ $`=`$ $`\left(T_c+V_{NN}\right)_{{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}\text{He}}\left(T_c+V_{NN}\right)_{\text{4}\text{He}}`$ (9)
$``$ $`4.7\text{MeV}.`$
This difference is considerably large despite the fact that the rms radius of $`N`$ from the CM of 4He for the $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He hardly changes from that for 4He. (Both radii are 1.5 fm.) Most of the change is due to a reduction of the energy expectation value of the tensor $`NN`$ interaction,

$$\left(V_{NN}(\text{tensor})\right)_{{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}\text{He}}\left(V_{NN}(\text{tensor})\right)_{\text{4}\text{He}}2.9\text{MeV}.$$

On the other hand, the tensor $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ transition part has a surprisingly large energy expectation value (about $`20`$ MeV). This large coupling energy makes $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He bound in spite of both the energy loss of $`\mathrm{\Delta }E_c`$ and the extremely high energy of the $`\mathrm{\Sigma }`$-component ($`\frac{H_\mathrm{\Sigma }}{P_\mathrm{\Sigma }}600`$ MeV).

The calculated wave function is divided into orthogonal components according to the total orbital angular momentum ($`L`$), the total spin ($`S`$), the core nucleus spin ($`S_c`$) and the core nucleus isospin ($`T_c`$). Table 4 displays the probability of each component for $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He. The table also lists the probability of $`S`$-state or of $`D`$-state for 4He. The sizable amount of probability of the $`\mathrm{\Sigma }`$-component is found in the $`D`$-state while the sum of $`D`$-state probabilities in the $`\mathrm{\Lambda }`$-component is slightly smaller than that for 4He. Moreover, though the presence of a $`\mathrm{\Lambda }`$ in 4He with the strong tensor $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ transition potential influences the structure of the $`D`$-state component and reduces the energy expectation value of the tensor $`NN`$ interaction, the large coupling energy $`V_{\mathrm{\Lambda }\mathrm{\Sigma }}`$ of the tensor part bears the bound state of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He instead.

In summary, we have made a systematic study of all $`s`$-shell hypernuclei based on ab initio calculations using $`YN`$ interactions with an explicit $`\mathrm{\Sigma }`$ admixture. The bound-state solution of $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He was obtained. As the Ref. Nogga  claimed, though there is none of the interaction models to describe very precisely the experimental $`B_\mathrm{\Lambda }`$’s, the five-body calculation convinced us that the anomalous binding problem would be resolved by taking account of the explicit $`\mathrm{\Sigma }`$ admixture. The $`B_\mathrm{\Lambda }`$ values, obtained by using SC97e(S), are the closest to the experimental values, among the $`YN`$ interactions employed in this study. A sizable amount of $`P_\mathrm{\Sigma }`$ was obtained, even for the $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He, in spite of the large excitation energy of 4He. The contribution of the energy from the tensor $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ coupling is quite large, and this coupling is considerably important to make $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He bound. This is novel finding in contrast with the Brueckner-Hartree-Fock calculation Akaishi . The present study for $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He is a first step toward a detailed description of light strange nuclear systems. The core nucleus, 4He, is no longer rigid in interacting with a $`\mathrm{\Lambda }`$ particle. A similar situation can occur for strangeness $`S=2`$ systems. Investigations into the strength of the $`\mathrm{\Lambda }\mathrm{\Lambda }`$ interaction based on the experimental data of the binding energy for double $`\mathrm{\Lambda }`$ hypernuclei (e.g. $`{}_{}{}^{6}{}_{\mathrm{\Lambda }\mathrm{\Lambda }}{}^{}`$He Nagara ) should take account of the energy reduction of the core nucleus ($`\mathrm{\Delta }E_c`$ for $`{}_{}{}^{6}{}_{\mathrm{\Lambda }\mathrm{\Lambda }}{}^{}`$He is expected to be larger than the present $`\mathrm{\Delta }E_c`$ for $`{}_{}{}^{5}{}_{\mathrm{\Lambda }}{}^{}`$He).

###### Acknowledgements.

We are thankful to Y. Fujiwara, K. Miyagawa, H. Kamada, K. Varga and S. Shinmura for useful communications. One of the authors (H.N.) would like to thank for JSPS Research Fellowships for Young Scientists. The calculations of $`A=4,5`$ systems were made using the RCNP’s SX-5 computer and KEK’s SR8000 computer.


format: text

text:
# Thermodynamic driving force of formation of coherent three-dimensional islands in Stranski-Krastanov growth

## Abstract

The formation of coherent three-dimensional islands in highly mismatched epitaxy is discussed in terms of the traditional concept of wetting. It is shown that the wetting layer and the 3D islands represent different phases which cannot be in equilibrium with each other. The transfer of matter from the stable wetting layer to the 3D islands is thermodynamically unfavored. The experimentally observed critical misfit for coherent 3D islanding to occur and the coexistence of pyramids with discrete heights of two, three, four… monolayers can be explained assuming that the thermodynamic driving force of formation of coherent 3D islands on the surface of the wetting layer of the same material is the reduced average adhesion of the islands to that layer and that the islands height is a discrete variable.

The growth of thin epitaxial films usually takes place far from equilibrium. Nevertheless, thermodynamic considerations are a necessary step for understanding of the process. Of particular interest is the thermodynamic driving force (TDF) which is responsible for one or another mechanism of growth. While this question is well understood in terms of wetting of the substrate by the overgrowth in the cases of island or Volmer-Weber (VW) growth and layer-by-layer or Frank-van der Merwe (FM) growth,Bauer ; Kern ; Gilmer  the Stranski-Krastanov (SK) growth (3D islands on top of a thin wetting layer) is far from being clarified. The reason is that the SK growth is in fact a growth of material A on the same material A, which thermodynamically requires the formation and growth of 2D rather than 3D islands. This is particularly true in the case of the coherent Stranski-Krastanov growth,Eagle  where dislocation-free 3D islands are strained to the same degree as the wetting layer.Eagle ; Mo ; Moison  This is the reason why it is widely accepted that the energy of the interfacial boundary between the 3D islands and the wetting layer is equal to zero.J1  Although this energy is expected to be small compared with that of the free crystal faces,Jacques  it should not be neglected since this is equivalent to the assumption that the islands wet completely the wetting layer. The latter rules out 3D islanding from a thermodynamic point of view.Sir

The need for a thermodynamic analysis arises also from the experimental observations of a critical misfit for coherent 3D islanding to occur,Xie ; Leo ; Walther ; Pinc  and the simultaneous presence of islands of different thickness which vary by one monolayer.Rud ; Col  The existence of a critical misfit, as well as of stable two, three or four monolayers thick islands, do not follow from the tradeoffJ1 ; Jacques

$`\mathrm{\Delta }EC^{}\gamma V^{2/3}C^{\prime \prime }\epsilon _0^2V`$ (1)
between the cost of the additional surface energy and the gain of energy due to the elastic relaxation of the 3D islands relative to the wetting layer ($`V`$, $`\gamma `$ and $`\epsilon _0`$ are the islands volume, the specific surface energy and the lattice misfit, respectively, and $`C^{}`$ and $`C^{\prime \prime }`$ are constants).

It was recently suggested that the TDF for coherent 3D islanding is the incomplete wetting of the substrate by the islands,Kor  rather than the elastic relaxation of the material in the islands. The incomplete wetting is due to the displacements of the atoms near the island edges from the bottoms of the corresponding potential troughs provided by the wetting layer. This results in a series of critical volumes at which the monolayer high islands become unstable against the bilayer islands, the bilayer islands against the trilayer islands, etc. The misfit dependence of the first critical size $`N_{12}`$ for the mono-bilayer transformation displays a critical behavior in the sense that coherent 3D islands can be formed at a misfit higher than some critical value. Below this value the film should grow in a layer-like mode until misfit dislocations are introduced to relieve the strain. However, the approximation used by the authors, which is based on the 1D model of Frenkel and Kontorova,Ratsch ; Frenkel  was unable to describe correctly the individual behavior of atoms inside each layer, since it assumes a potential with a period given by the average of the separations of atoms (considered frozen) in the layer underneath. Although this model gives qualitatively reasonable results concerning the energy of the islands, it is inadequate to calculate, in particular, the average adhesion energy of the islands to the wetting layer.

In the present report we recollect some simple thermodynamic aspects of the epitaxial morphology based on the traditional concept of wetting and consider the coherent SK growth from this point of view. The same concepts were in fact advanced by Stranski in his model, admittedly very peculiar, of a monovalent ionic crystal K+A- on the surface of an isomorphous bivalent crystal K2+A2-.Stran1  We then support our thermodynamic considerations by numerical calculations making use of a simple minimization procedure on the same atomistic model in $`1+1`$ dimensions (length + height) as in Ref.(Kor, ). The 3D islands are represented by linear chains of atoms stacked one upon the other,Ratsch  the islands height being thus considered as a discrete variable which increases by unity from one. The latter is of crucial importance as the aspect ratio of the 3D islands is usually of the order of 0.1 and the height is of the order of 10 monolayers.Mo ; Moison  The atoms interact through an anharmonic Morse potential

$$V(x)=V_o[e^{12(xb)}2e^{6(xb)}].$$
(2)
The total interaction energy as well as its derivatives with respect to the atomic coordinates, i.e. the forces, are calculated. Relaxation is then performed iteratively by allowing the atoms to displace in the direction of the forces until these fall below some negligible cutoff value. We consider interactions in the first coordination sphere in order to mimic the directional bonds that are characteristic for semiconductors.Tersoff  Inclusion of further coordination spheres alters only minimally the numerical results. The substrate (the wetting layer) is assumed to be rigid.

A mother phase (a vapor) and a new phase (e.g. a strained planar film or unstrained 3D crystals) are in equilibrium with each other when their chemical potentials are equal. A transition from one phase (mother or new) to another takes place when the chemical potential of one of the phases becomes smaller than that of the other. The TDF for this transition is the difference of the chemical potentials of both phases at the given pressure and temperature. The TDF which determines the occurrence of one or another mechanism of epitaxial growth (growth from vapor of a strained 2D layer or 3D islands) is the difference $`\mathrm{\Delta }\mu =\mu (n)\mu _{3\mathrm{D}}^0`$ of the chemical potential, $`\mu (n)`$, of the overlayer which depends on the film thickness measured in number $`n`$ of monolayers counted from the interface, and the chemical potential, $`\mu _{3\mathrm{D}}^0`$, of the bulk 3D crystal of the same material.Kern ; Gilmer  The thickness dependence of $`\mu (n)`$ originates from the thickness distribution of the misfit strain and, on the other hand, from the interaction between the deposit and the substrate, which rapidly decreases with the distance from the interface ($`E_{\mathrm{AB}}E_{\mathrm{AA}}`$).Kern ; Mark1

If we deposit a crystal A on the surface of a crystal B, $`\mathrm{\Delta }\mu `$ can be written in terms of the interatomic energies per atom, $`E_{\mathrm{AA}}`$ and $`E_{\mathrm{AB}}`$, required to disjoin a half-crystal A from a like half-crystal A and from an unlike half-crystal B, respectively:Kaisch

$$\mu (n)=\mu _{3\mathrm{D}}^0+[E_{\mathrm{AA}}E_{\mathrm{AB}}(n)]=\mu _{3\mathrm{D}}^0+E_{\mathrm{AA}}\Phi .$$
(3)
The adhesion energy $`E_{\mathrm{AB}}`$ includes in itself the thickness distribution of the strain energy and the attenuation of the bonding with the substrate. $`\Phi =1E_{\mathrm{AB}}/E_{\mathrm{AA}}`$ is the adhesion parameter which accounts for the wetting of the substrate by the overgrowth. Eq. (3) is equivalent to the familiar 3-$`\sigma `$ criterion of Bauer.Bauer ; Mark1

As follows from (3) the parameter $`\Phi =\mathrm{\Delta }\mu /E_{\mathrm{AA}}`$ is equal to the TDF for occurrence of one or another mode of growth relative to the cohesive energy $`E_{\mathrm{AA}}`$. In the two limiting cases of VW ($`0<\Phi <1`$) and FM growth ($`\Phi 0`$, $`\epsilon _00`$), $`\mathrm{\Delta }\mu `$ tends asymptotically with increasing film thickness to zero from above and from below, respectively, but changes its sign in the case of SK growth ($`\Phi <0`$, $`\epsilon _00`$), as shown in Fig. 1.Gilmer ; Mark1  Consider now Fig. 1 in terms of equilibrium vapor pressures instead of chemical potentials. Although the connection is straightforward, ($`\mu \mathrm{ln}P`$), such a consideration gives a deeper insight into the problem.Stran1  Thus, as long as $`\mu (n)<\mu _{3\mathrm{D}}^0`$ a thin planar film can be deposited at a vapor pressure $`P`$ that is smaller than the equilibrium vapor pressure, $`P_0`$, of the bulk crystal, but is larger than the equilibrium vapor pressure $`P_1`$ of the first monolayer, i.e. $`P_1<P<P_0`$. In other words, a planar film can be deposited at undersaturation $`\mathrm{\Delta }\mu =kT\mathrm{ln}(P/P_0)`$ with respect to the bulk crystal. The formation of 3D islands ($`\mu (n)>\mu _{3\mathrm{D}}^0`$) requires $`P>P_0`$, or a supersaturation with respect to the bulk crystal.

Applying the above considerations to the SK growth leads unavoidably to the conclusion that the 3D islands and the wetting layer represent necessarily different phases and thus have different chemical potentials. The reason is that the two phases are in equilibrium with the mother phase (the vapor) under different conditions which never overlap. The wetting layer can be in equilibrium only with an undersaturated vapor phase ($`P<P_0`$), while small 3D islands can be in equilibrium only with a supersaturated vapor phase ($`P>P_0`$). The dividing line is thus $`\mathrm{\Delta }\mu =kT\mathrm{ln}(P/P_0)=0`$ at which the wetting layer cannot grow thicker and the 3D islands cannot nucleate and grow. Hence the wetting layer and the 3D islands can never be in equilibrium with each other.

It follows from the above that the derivative, $`d\mathrm{\Delta }E/dV`$, of the energy of the 3D islands relative to that of the wetting layer, gives the difference of the chemical potentials of the wetting layer (the dashed line in Fig. 1) and the chemical potential of the 3D islands. In other words, it represents the difference of the supersaturations of the vapor phase with respect to the wetting layer and the 3D islands. As the thickness and the energy of the wetting layer depend on the misfit it would be more suitable if one chooses as a reference the bulk crystal rather than the wetting layer.Muller  Transfer of material from the stable wetting layer ($`\mu _{\mathrm{WL}}<\mu _{3\mathrm{D}}^0`$) to the 3D islands is connected with increase of the free energy of the system, and therefore, is thermodynamically unfavored. A planar film thicker than the stable wetting layer is unstable and the excess of the material can be transferred to the 3D islands if the necessary thermal activation exists.

We focus our attention on the adhesion between the 3D islands and the wetting layer. Fig. 2 is an illustration of the difference (and resemblance) between the classical and the coherent SK mode. In the calculations, the sizes of the base chain of the island (30 and 34 atoms, respectively) have been chosen just below and above the critical size for introduction of misfit dislocations at the given misfit of 7%. As seen, in both cases the 3D islands loose contact with the wetting layer. The vertical displacements are largest at the chain’s ends in the coherent SK mode and around the dislocation cores in the classical case, but the physics is essentially the same. The mean adhesion parameter $`\Phi `$ increases with the islands’ height and saturates beyond several monolayers (Fig. 3). In our model, $`\Phi `$ is calculated as the adhesion energy between island and wetting layer at the given misfit (7%) minus the corresponding value for zero misfit. It can be seen that the compressed overlayers exhibit a greater tendency to coherent SK growth than expanded ones as expected, due to the anharmonicity of the potential (2).

The same tendency is seen more clearly in Fig. 4 which shows in fact the dependence of the TDF for formation of coherent 3D islands, $`\mathrm{\Delta }\mu `$, on the lattice misfit. The latter remains close to zero for expanded overlayers but increases steeply beyond approximately 5% in compressed overlayers. This behavior agrees well with the misfit dependence of the critical size, $`N_{12}`$, for the mono-bilayer transformation to occur, as shown in Fig. 5, where a steep rise of $`N_{12}`$ with decreasing absolute value of the misfit is observed only in compressed overlayers. Note that it is less sharp for expanded ones opposite to earlier finding.Kor

We discuss now the discrete character of the height of the 3D islands. The experimentally observed volume of the quantum dots varies roughly from 20000 to 50000 atoms.Moison ; Leo  Typical values of the aspect ratio of the islands height and half-base are of the order of 0.1.Moison ; Voigt  Thus a pyramid with a base edge of 100 atoms and aspect ratio 0.1, and containing 22000 atoms, is only 5 monolayers high. The addition of 14400 atoms (a new base plane of $`120\times 120`$ atoms) requires only one more atomic plane. Calculations of the energy of islands having a shape of a frustum of a pyramid are usually performed assuming implicitly that the lengths of the lower, $`R`$, and the upper, $`R^{}`$, bases, and in particular the height $`h`$, are continuous variables. Eq. (1) is obtained by using the Tersoff approximation which neglects the gradient of strain in a direction normal to the surface plane together with $`hR`$,J1  and assuming $`hc`$ where $`c`$ is the atomic distance.Jacques  This would be correct if the crystals contain at least several million of atoms.

We conclude that the wetting layer and the 3D islands represent different phases which cannot be in equilibrium with each other, and the SK morphology is a result of the replacement of one first order phase transition (vapor – wetting layer) by another first order transition (vapor – 3D islands). The transfer of matter from the stable wetting layer to the 3D islands is thermodynamically unfavored. The experimental observations mentioned above can be explained on the base of two assumptions: the thermodynamic driving force for the coherent 3D islanding is the incomplete wetting and the height of the 3D islands is a discrete variable varying by one monolayer. This leads to the results that (i) monolayer high islands with a critical size appear as necessary precursors for 3D islands, (ii) the 2D-3D transition takes place through a series of intermediate states with discretely increasing thickness that are stable in separate intervals of volume (see Refs. (Rud, ; Col, )) (iii) there exists a critical misfit below which coherent 3D islands are thermodynamically unfavored (see Refs. (Xie, ; Leo, ; Walther, ; Pinc, )) and the misfit is accommodated by misfit dislocations at a later stage of the growth. Compressed overlayers show a greater tendency to 3D clustering than expanded ones, in agreement with experimental results.Xie  Result (i) explains readily why the volume distribution of InAs/GaAs self-assembled quantum dots agrees well with the scaling functions for two-dimensional submonolayer homoepitaxy model.Ebiko

The authors are indebted to the Instituto Universitario de Ciencia de Materiales “Nicolás Cabrera” for granting research visits which enabled scientific collaboration. This work was supported by the Spanish CICyT through project Nr. MAT98-0965-C04-02.


text:
# HST Survey of Clusters in Nearby Galaxies. II. Statistical Analysis of Cluster Populations

## 1 Introduction

As with stars, the star clusters of a galaxy provide direct fossil evidence of previous star formation events. Compact clusters span the complete range of ages in the universe, from very recent star formation (30 Doradus) to the oldest components of galaxies (globular clusters). Integrated photometry of a single cluster can provide an age and initial mass of that cluster, thus allowing one to trace the cluster formation history in a similar way that blue helium-burning stars trace the star formation history. However, the effects of stellar evolution, cluster lifetimes, differential reddening, etc. prevent such a simple determination of physical parameters such as the cluster initial mass function and cluster formation rate.

Because of the difficulties in determining the physical parameters, work on cluster populations have used relatively crude statistical comparisons of observed properties. The simplest is a measurement of the combined light of detected clusters, which can be compared with the brightness of the parent galaxy to determine a cluster formation rate parameter (Larsen & Richtler, 2000). The advantage of such a technique is its relatively easy implementation and minimal amount of required data. For example, WLM (a Local Group) galaxy has only one globular cluster. One would not be able to determine a luminosity function, star formation history, or any other detailed parameters from this object; however a cluster frequency can be determined very accurately.

A more sophisticated technique involves an examination of the cluster luminosity function. Whitmore et al. (1999) have made very effective use of this in their study of the Antennae, where there are many thousands of clusters to build up accurate statistics. While a luminosity function contains considerably more information than does a cluster frequency, the complications listed above make the physical interpretation of the observed luminosity function extremely difficult, except under special conditions such as all clusters having the same age.

What is largely lacking is a comprehensive statistical technique to interpret the observational data in terms of physical properties. Much progress has been made in the ability to determine star formation histories from stellar color-magnitude diagrams Dolphin (1997, 2002); little progress has been made in the equivalent study of cluster populations.

This paper is the second in our series on a survey of the cluster populations of nearby galaxies. Dolphin & Kennicutt (2002; hereafter Paper I) demonstrated our procedure for the identification and photometry of clusters in our data, which consist of WFPC2 observations of galaxies between $`4`$ and $`25`$ Mpc. In this paper, we are examining statistical methods of cluster populations with which we can constrain physical properties.

In section 2, we examine the many factors that affect the observed cluster population of a galaxy and describe a technique that will incorporate all of these to produce a synthetic cluster population for any set of parameters. In section 3, we show how synthetic cluster populations and real cluster populations can be statistically compared to maximize the physical information; this technique is tested using simulated data in section 4. Finally, we will conclude our study of NGC 3627 (begun in Paper I) with a demonstration of our statistical technique in Section 5.

## 2 Synthetic Cluster Populations

The key to accurate analysis of a galaxy’s cluster population is an understanding of how various physical characteristics (such as cluster size) transform into observational quantities (such as sharpness). Most observational quantities are functions of more than one physical characteristic. For example a cluster’s color is related to its age, metallicity, extinction, and the random selection of stellar masses; thus a simplistic attempt to measure one of those quantities using color is likely to fail. In the same way, completeness is a function of color, magnitude, and cluster radius; an attempt to use the observed luminosity function to determine the cluster mass function is likewise inaccurate unless color and radius distributions are taken into account. We will examine each of the physical parameters in this section. Some will be seen to be very important; others will be less so.

Our analysis approach will center on generating and analyzing synthetic cluster populations. Because our photometry provides three variables (magnitude, color, and sharpness), we will refer to color-magnitude-radius diagrams (CMRDs) instead of color-magnitude diagrams (CMDs). In this section we describe the ingredients of a synthetic CMRD. Our approach will be parallel to the star formation history measurement technique laid out by Dolphin (1997, 2002); we will generate “noiseless” CMRDs that account for all possible combinations of age, metallicity, mass, and photometric effects in order to make statistical comparisons between these and the observed data.

### 2.1 Cluster Formation History

The primary influences on a galaxy’s cluster population are its star formation history (mass of stars formed as a function of age and metallicity), cluster formation efficiency (mass of clusters formed divided by mass of stars formed), and chemical enrichment history. However, constraints for these quantities are very poorly-determined. The present-day mean metallicity is $`[\text{M/H}]=+0.2\pm 0.2`$ (Kennicutt et al., in prep) from nebular spectroscopy. We can reasonably assume that the mean metallicity has slowly and consistently increased from metal-poor globular abundances (for which we use $`[\text{M/H}]=2.3`$, the most metal-poor of the Girardi et al. (2000) isochrones) at ancient ages to the present-day value, which allows us to estimate an age-metallicity relation. We note that the metallicities at older ages are not as important, since clusters dim so rapidly that most of what we observe is younger than $`2040`$ Myr. Since there is some scatter in the metallicities of HII regions within a galaxy, we must also allow for scatter around our assumed age-metallicity relation. We adopt an rms scatter of $`\pm 0.13`$ dex, based on the Galactic disk study of Rocha-Pinto et al. (2000).

The star formation rate and cluster formation efficiency, however, are nearly completely unconstrained from the $`V`$ and $`I`$ imaging alone. However, the timescale over which clusters are visible is extremely brief. Our NGC 3627 data (Paper I, figure 15) shows that even a cluster with initial mass of $`10^5M_{}`$ would fade below the bottom of our observed CMRD in $`2040`$ Myr. We therefore will make the assumption that the cluster formation rate (the product of the star formation rate and cluster formation efficiency) has been constant over the observable period. This assumption would only be invalid in the case where a short burst created a large number of clusters in the same evolutionary phase; such a situation would be obvious from an examination of the CMRD, however, as all clusters formed in a burst will have roughly the same color. We recognize that, during the epoch of globular cluster formation, the cluster formation rate was likely much different than what is observed presently. Thus our synthetic cluster CMRD will consist of two epochs of cluster formation, with a dividing line of 10 Gyr. These two populations will be distinct in broadband $`VI`$ photometry, as the old population consists entirely of red ($`VI>0.8`$) clusters while the younger population contains a mixture of colors.

In the case of NGC 3627 (analyzed below in Section 5), we indeed find that the form of our cluster formation rate is unimportant, and that the 2-parameter star formation rate fits the data as well as a more sophisticated function. We also created synthetic CMRDs using a variety of metallicity enrichment laws (with the same current metallicity), all of which produced similarly-good fits to the data.

We therefore characterize the cluster formation history with four parameters: cluster formation rates at young and ancient ages, an age-metallicity relation, and scatter around the age-metallicity relation. However, while all of these ingredients are necessary to produce a synthetic CMRD, the only parameter that significantly affects the CMRD is the recent cluster formation rate.

### 2.2 Initial Mass and Radius Distributions

Once we have determined (or assumed) the cluster formation rates and the age-metallicity relation, we know the mass of clusters formed at each age and metallicity. To produce a synthetic CMRD, however, we need the actual number of clusters formed.

We first divide the mass of clusters at one age and metallicity according to an assumed cluster initial mass function. The function can be anything; we will assume a power law of the form

$$\frac{d\mathrm{log}(dN/d\mathrm{log}M)}{d\mathrm{log}M}=\alpha ,$$
(1)
or

$$\frac{dN}{d\mathrm{log}M}M^\alpha .$$
(2)
As this function diverges at low mass for $`\alpha 1`$, we arbitrarily choose a low-mass cutoff of 1 solar mass for now, and we hope that we can constrain the cluster initial mass function at low mass from the closer galaxies in our sample.

The second division of clusters will be according to cluster radius. This is extremely difficult, as our observed data have very strong selection effects in radius: small clusters will be classified as stars while large clusters will be too difficult to photometer or may be identified as small associations and thus rejected. Because we are sensitive to only a limited range of core radii ($`0.24`$ pc in our NGC 3627 data), we will initially assume a flat distribution of cluster radii,

$$\frac{dN}{dr}=constant.$$
(3)
We will be able to constrain this further by comparing the sharpness distributions of the observed and synthetic CMRDs.

With the clusters now divided by age, metallicity, initial cluster mass, and core radius, one other factor affects the physical properties of the cluster population: lifetimes. However, in terms of generating a synthetic CMRD, we consider this effect unimportant due to the very short timescales over which clusters are visible above our photometric cutoffs. With our limiting magnitude of $`V=24`$, we only require that a cluster with an initial mass of $`10^4M_{}`$ lives $`1015`$ Myr and that one with a mass of $`10^5M_{}`$ lives $`3040`$ Myr. (That the upper end of the mass function is unaffected by disruption is demonstrated by the analytical models of Fall & Zhang 2001.) As cluster lifetimes are expected to be at least an order of magnitude beyond these values (Kroupa, Aarseth, & Hurley, 2001; Fall & Zhang, 2001), we can thus ignore the effect of cluster disruption.

We have thus added two additional parameters to those listed at the end of the previous section: the cluster initial mass function and radius distribution.

### 2.3 Intrinsic Observed Properties

Accounting for the factors listed in the previous two sections, we now know the physical properties (age, metallicity, mass, and core radius) of the cluster system. These must now be converted into the intrinsic observable quantities of magnitude, color, and angular size. The final quantity is trivially determined from the core radius and distance; we address the others in this section.

The magnitude and color of a cluster are generally determined using stellar evolution models (Bruzual & Charlot, 1993; Girardi et al., 2000). We will adopt that approach here as well, using the models of Girardi et al. (2000) and the revised stellar IMF of Kroupa (2001) to determine integrated magnitudes and colors. We acknowledge that the theoretical stellar models are not perfect; however we believe that they are sufficiently accurate that integrated magnitudes of clusters should be reasonably correct, even if a few evolutionary phases are modeled inaccurately. At any rate, as no empirical database of integrated cluster magnitudes exists spanning the entire range of ages, metallicities, and filter passbands that we will need, theoretical models must be used.

There is an additional complexity that we must account for. Because the division of a certain amount of mass into stars does not always produce the same set of star masses, there is scatter inherent in the cluster colors. In other words, clusters of the same age and metallicity do not have exactly the same color and magnitude. We therefore used Monte Carlo simulations to determine the color and magnitude distribution of clusters with a wide range of age, metallicity, and initial mass; the full distribution is added to our CMRD rather than merely the mean magnitudes and colors.

For example, the mean cluster of solar metallicity, age 10 Myr, and initial mass of $`10^5M_{}`$ has magnitudes of $`M_V=11.89`$ and $`VI=0.71`$. Because of noise in the distribution of stars, however, this population of clusters appears as a “smear” with $`\sigma _V=0.08`$ magnitudes and $`\sigma _{VI}=0.07`$ magnitudes. This factor becomes extremely important for young clusters with lower masses, in which case the scatter from the distribution of stars becomes much larger than the photometric errors.

Once the distribution of magnitudes and colors is determined, we adjust the magnitudes for distance and global extinction. For NGC 3627, we determined the distance from Cepheids and adopted the foreground extinction of $`A_V=0.11`$, calculated from the maps of Schlegel, Finkbeiner, & Davis (1998). Because most (if not all) the galaxies in our sample have ongoing star formation and thus dust, we must also account for differential extinction. We therefore apply a range of extinction values, requiring the lowest value to equal the foreground extinction and the mean value to equal the extinction measured from the Cepheids. We will begin by assuming a flat distribution of extinction values; this can be constrained slightly by the data.

### 2.4 Observational Errors and Selection Effects

The final ingredient needed to generate a synthetic CMRD is the artificial cluster tests. These were described in full in Paper I; we will summarize here. Our artificial cluster tests are identical in concept to artificial star tests: we add artificial clusters to the image (using an assumed cluster profile convolved with the stellar PSF), photometer them, and apply the same magnitude, sharpness, and roundness cuts that were applied to the real data. This produces an artificial cluster library that characterizes not only the completeness fraction, but also the recovered magnitude, color, and sharpness as a function of input magnitude, color, and angular size.

Using this library, we can thus reproduce the photometry of the galaxy being studied, as well as exactly reproduce our selection effects. This is the reason why we chose very conservative magnitude and sharpness limits in Paper I; accounting for false positive detections is much more difficult than simulating selection effects.

### 2.5 Building a Synthetic CMRD

With all of the ingredients described in this section, we are able to generate a synthetic cluster CMRD. Summarizing the parameters listed above, the following pieces of information are all that is required:

1. Young and ancient cluster formation rates
2. Age-metallicity relation and scatter
3. Cluster mass and radius distributions
4. Distance, extinction, and differential extinction
5. Synthetic models from isochrones, including intrinsic scatter in colors and magnitudes
6. Artificial cluster library

Before we begin an analysis of a galaxy, the final two items are already calculated. The synthetic models are pre-computed, while the artificial cluster library is part of our photometry process. The distance is measured from Cepheid photometry, as is an initial guess of the extinction. We assume an age-metallicity relation, and assume functional forms for the cluster initial mass function, radius distribution, and differential extinction. An initial guess as to the mass function can be taken from the luminosity function. As noted above, the resulting CMRD is almost completely insensitive to the age-metallicity relation, provided that the current metallicity is fixed to the observed value. This leaves only two completely unknown parameters (the young and ancient cluster formation rates) and three constrainable parameters with initial guesses (extinction, mass function, and radius distribution).

## 3 Cluster Population Statistics

The previous section describes how we generate a noiseless synthetic CMRD, which samples all possible ages, metallicities, masses, radii, extinction values, and photometric effects. As described in section 2.5, all but three parameters are known or assumed before we begin our analysis: the young cluster formation rate, the ancient ($`>10`$ Gyr) cluster formation rate, and the slope of the cluster initial mass function. These three values can be determined by a statistical comparison of the observed and synthetic CMRDs.

Computationally, we characterize each CMRD as a three-dimensional matrix, with the dimensions corresponding to magnitude, color, and sharpness. Within each element of the matrix is the number of observed or synthetic stars with that magnitude, color, and sharpness. Because of the nature of WFPC2 data, various parts of the image will have different statistical characteristics. For NGC 3627, for example, WFC2 contained severe crowding (because of its proximity to the center of the galaxy) not seen in the other chips. Additionally, the PC always has different characteristics than the other chips because of its different plate scale. Thus we produce multiple CMRDs, one corresponding to each region. For most galaxies, we will use one CMRD per chip. When studying galaxies with photometry in three or more bands, we will again multiply the number of CMRDs to allow the simultaneous examination of, for example, a (B, B-V) CMRD and a (V, V-I) CMRD.

We judge the level of agreement between the observed and synthetic CMRDs by the probability that the observed data were drawn from the synthetic CMRD. Because our synthetic CMRD is noiseless (rather than randomly populated), we can use the Poisson probability law:

$$P=\underset{i}{}\frac{e^{m_i}m_i^{n_i}}{n_i!},$$
(4)
where $`P`$ is the probability, $`i`$ represents the CMRD element, $`m`$ is the synthetic (model) CMRD, and $`n`$ is the observed CMRD. Maximizing this value for an unchanging observed CMRD is the equivalent of minimizing

$$2\underset{i}{}m_in_i+n_i\mathrm{ln}\frac{n_i}{m_i},$$
(5)
the Poisson equivalent of $`\chi ^2`$. By determining the combination of cluster formation rates and mass function slope that minimizes this equation, as well as the acceptable range of values surrounding the minimum, we can thus determine these values empirically. We also anticipate being able to constrain the cluster radius distribution. Once we have compiled sufficient data, we will also attempt to constrain the lower-mass cutoff (or any slope change) in the mass function.

Details of the statistics, including a derivation of this formula, an explanation of why $`\chi ^2`$ cannot be used with Poisson-distributed data, importance of using noiseless synthetic CMRDs, and techniques for evaluating the quality of the best fit and uncertainties in the parameters are found in Dolphin (2002); we do not wish to duplicate the lengthy discussion of statistics here and instead refer the reader to that paper.

## 4 Synthetic Galaxy

Before concluding our example analysis of NGC 3627, we wish to demonstrate the effectiveness of our analysis technique via application to a simulated galaxy. To make the results of this exercise directly comparable to our NGC 3627 study, we will use the parameters producing the best fit to NGC 3627’s cluster population. These values are shown as the “input” values in Table 1.

### 4.1 Fit Accuracy

Our first test was to determine the parameters that produced the best fit to the synthetic data. A comparison between these values and those used to create the synthetic data will test our technique for any biases. While some difference is expected (because the simulated observed CMRD was randomly populated), one expects that the differences between the recovered and input values will be less than the uncertainties. Indeed, this was true, as is seen in Table 1. Additionally, the recovered cluster formation rate differed from the input rate by less than 1%. Thus, with all values being measured within their uncertainties, we have no reason to suspect biases in the method.

### 4.2 Correlated Parameters

Our second test is a search for correlated errors in the parameters, which was done by forcing each of the parameters, in turn, to a “bad” value ($`1\sigma `$ from the best value) and looking for any effects on the best fits of the other parameters. Our expectation is that the mean value and width of the radius distribution will be correlated with each other (since both are primarily related to the cluster sharpness), but that the cluster IMF slope (primarily a function of $`V`$ magnitude) and extinction distribution (primarily a function of $`VI`$ color) will be independent. Table 2 shows the results of this test, with deviations given in units of $`\sigma `$. One can easily discern the correlations between the two radius parameters, as well as their independence from the other two parameters. What is more surprising is the apparent correlation between the cluster IMF slope and the extinction distribution. Since no such correlation should exist (a pure power-law luminosity function will keep the same slope, regardless of extinction), we assume this to be a numerical instability created by the poorly-constrained nature of the extinction, rather than a true correlation. We note that the change in the mass function slope created by the extinction change is extremely small ($`0.05`$); it appears as a large number in the table because the uncertainty in the measured value is $`\pm 0.07`$.

## 5 NGC 3627

With the mechanism now in place for analyzing the cluster population, we return to the study of NGC 3627 (M66) begun in Paper I. As noted in Paper I, this galaxy is part of an interacting system. We found 528 clusters, which account for $`1`$% of the galaxy’s luminosity. By re-reducing the data with a newer version (1.1) of HSTphot (Dolphin, 2000), the number of clusters is now 506. We have also conducted extensive artificial cluster tests, in which clusters of known integrated magnitudes and radii are inserted and photometered; this allows us to accurately model observational effects such as incompleteness, photon noise, crowding, and photometric biases.

For our analysis, we have reduced the data to two sets. First is the set of 506 observed clusters, for which we need only the observed $`V`$ and $`I`$ magnitudes, as well as the sharpness measurement. Photometry uncertainties are not needed, as these are simulated using the artificial clusters. The second data set is the artificial cluster library, which contains input integrated $`V`$ and $`I`$ magnitudes and core radii, as well as the observed properties (or a flag if the cluster was either not recovered or did not meet our cluster selection criteria). The CMRD modeling then proceeded according to the steps described in the previous sections.

### 5.1 Fit Quality

Although this modeling is not quite as sophisticated as the stellar population CMD-fitting method of Dolphin & Kennicutt (2002) due to a number of approximations in our synthetic CMRD construction, we believe that it is useful to look at the fit quality to determine how well we matched the observed data. The maximum likelihood fit parameter (equation 5) of our best solution was 924.7. For comparison, the typical data set randomly drawn from our best model would produce a fit parameter of 850.9, with a $`1\sigma `$ scatter of 38.2. Thus the match between the observed data and our best model is roughly 2 $`\sigma `$ from ideal. Given the approximations involved in using functional forms for differential extinction and the mass and radius distributions, we consider this a good fit.

A comparison between the top and bottom panels of Figure 1 shows that our best model matched the observed sharpness distribution very well. A comparison of the CMDs is given in Figure 2. The top panels show the observed (left) and model (right) CMDs in greyscale; the break at V=23.5 is caused by the higher cutoff magnitude in WFC2 that was imposed to combat contamination. (Since the magnitude limits were also applied to the synthetic data, it is possible to include photometry of varying depth without any extra effort.) The bottom left panel shows the difference of the top panels, with darker squares signifying more observed than synthetic clusters and lighter squares the reverse. The important panel of Figure 2 is the bottom right panel, which shows the significance of the differences (essentially a $`\chi `$ measurement, though using Poisson statistics). In no part of the CMD is that fit worse than $`3\sigma `$, and only two CMD regions are worse than $`2\sigma `$. (Note that the CMRD as fit by our algorithm is truly a three-dimensional diagram; we have divided it into a CMD and a radius histogram in this section merely for ease of display.)

### 5.2 Cluster Initial Masses

Of the free parameters that enter into the solution, two are quite well-constrained. The first is the cluster initial mass function slope. Because the luminosity function of a population of coeval clusters is nearly the same as its mass function, and because the sum of many power laws with identical slopes is another power law with that slope, the luminosity function of a cluster population is a power law with the same slope as its mass function. Two factors complicate this slightly. First, the random division of a cluster’s mass into stars will cause the mass-to-light ratios of coeval clusters to be slightly different. As the scatter is larger for low-mass clusters, and such clusters are more common, the luminosity function is systematically flattened slightly. The second complication comes from incompleteness, which is a function of magnitude, color, and cluster size. Since completeness is lower at fainter magnitudes, this again flattens the luminosity function. Both of these are second-order effects, however; the relationship between the mass function and luminosity function is very strong and thus it is the best-constrained of the free parameters.

Using a trial-and-error approach (creating a variety of model CMRDs and determining which fits the best), we measure a mass function slope of $`\alpha =1.50\pm 0.07`$. This is slightly steeper than the observed luminosity function slopes of $`1.43\pm 0.12`$ (all clusters) and $`1.22\pm 0.23`$ (blue clusters) calculated in Paper I; this was expected since those values did not incorporate completeness corrections.

### 5.3 Cluster Radii

The second well-constrained input parameter is the distribution of cluster radii. As shown in Figure 3, the sharpness value measured by HSTphot is an extremely accurate measurement of the cluster’s radius. While there is some scatter in the diagram, the relation is tight in all chips, with sharpness proportional to the square root of the radius.

While the scatter in the sharpness vs. core radius relation and the strong influence of radius on completeness prevent one from determining the intrinsic radius distribution from the CMRD in the way that one can determine a luminosity function, it is not difficult to constrain the radius distribution by comparing the observed and synthetic CMRDs. Figure 1 shows the distribution of observed sharpness values, as well as those created by two different models. Because the middle panel (computed using a flat distribution of core radii) fits the observed data poorly, we can conclude that the shape of the observed distribution is not solely caused by selection effects and incompleteness. Instead, we are able to solve for the radius distribution.

The best-fitting sharpness distribution is shown in the bottom panel of Figure 1, and was created with a Lorentzian distribution of core radii centered at $`r_c=1.53\pm 0.15`$ pc and with a HWHM of 0.88 pc. The mean core radius is consistent with the results of Elson, Fall, & Freeman (1987), who studied the profiles of 10 young LMC clusters. The mean core radius of their sample of clusters was 2.19 parsecs, though the distribution is skewed by a few large clusters. The peak of their distribution is near a core radius of 1.7 parsecs with the values for most clusters falling between 1.3 and 2.

### 5.4 Extinction

Unlike the distributions of cluster mass and radii, the extinction profile is not readily apparent from the CMRD. We therefore took our initial guess of the extinction based on three observable values. The lower extinction limit was set, based on the maps of Schlegel, Finkbeiner, & Davis (1998), at $`A_V=0.11`$. The upper limit was determined by our reddest clusters; with no bright clusters redder than $`(VI)=1.8`$, there was no need to include extinction values greater than $`A_V=2.0`$. The final initial constraint was our observed mean Cepheid extinction values of $`0.58\pm 0.12`$ using the Madore & Freedman (1991) period-luminosity relation and $`0.76\pm 0.11`$ using the Udalski et al. (1999) relation.

Using these loose constraints, we attempted to fit the data using a variety of distributions. Because of the weak relationship between the extinction distribution and the shape of the CMRD, this solution was very poorly constrained. (The extinction determination was constrained by the color distribution of the clusters, but cluster formation history also affects color distribution and therefore the two factors are partially degenerate. Note that many of the galaxies we plan to study do have images in three broadband filters; in these cases the degeneracy can be more easily broken.) The function producing the CMRD best matching the data was a power law of the form

$$\frac{dN(A_V)}{dA_V}=A_V^{0.85\pm 0.36}.$$
(6)
Within our limiting extinction values, this distribution has a mean extinction of $`A_V=0.71\pm 0.15`$, which is in agreement with both measurements from the Cepheids.

### 5.5 Cluster Formation Rates

The final parameters in the solution are the young and ancient cluster formation rates. Unfortunately, the values we calculate are of limited value for a variety of reasons. Most importantly, from a single set of deep WFPC2 images, we have no reliable constraint on the mass function at high (due to WFPC2 saturation) or low masses (due to a loss of signal). Thus, while we were able to accurately measure the slope and scale of the mass function in the range where these data are sensitive, any estimate of the total mass contained in clusters is uncertain by many orders of magnitude. Note that calculating the combined cluster mass for the range of masses to which we are sensitive is not productive, since our galaxy sample includes objects spread over 5 magnitudes in distance modulus and WFPC2’s dynamic range (decreased by the high background under clusters) is only about 5 magnitudes. Therefore the set of clusters we can see at $`2530`$ Mpc do not overlap from those we can see at $`2.53`$ Mpc.

Our future papers will solve this problem in two ways. The direct approach is to observe galaxies at larger and smaller distances; examination of their luminosity functions will allow us to better-understand the cluster mass function at high and low mass. Even should that not provide adequate constraints, we can accurately determine relative cluster formation rates of a sample of galaxies by using the same mass function for all.

Secondly, our functional form of the extinction distribution is almost certainly an approximation. Since both the extinction and the cluster formation history affect the colors, any error in the measured extinction law will affect the derived cluster formation rates. Thus we do not feel that the cluster formation history of any one system can be accurately measured unless a third band of photometry is obtained.

## 6 Summary

We have presented a statistical system that can be used to study cluster populations. Our approach is to build synthetic cluster CMRDs for a variety of cluster population parameters, rather than attempting to take measurements directly from the observed data. In short, we generate a set of synthetic CMRDs based on a variety of cluster population properties and use a Poisson likelihood statistic to determine which is the best. Because of a variety of factors (completeness that varies as a function of color, magnitude, and radius; biases in the photometry; Poisson noise in the stellar composition of same-mass clusters; etc.), we believe that it is impossible to make definitive claims based on an examination of only the observed CMD; this statistical approach provides the capability of accounting for all such factors.

This approach is tested on the NGC 3627 data set from Paper I. We found that we were able to create a model that matched the observed data very well, both in reproducing the observed distribution of sharpness and the observed CMD.

The cluster mass function and radius distributions were very well constrained from the data. The cluster mass function is primarily determined by the cluster magnitudes, while the radius distribution is constrained by the sharpness measurements. Thus we can accurately determine a mass function slope of $`\alpha =1.50\pm 0.07`$ and a distribution of core radii centered on $`r_c=1.53\pm 0.15`$pc. Less certain were our measurements of the extinction distribution and cluster formation history, as both affect the colors. These variables will be more strongly constrainted in our future work, which will include galaxies observed through three broadband filters. Nevertheless, we were able to measure a mean extinction that was consistent with that measured from the Cepheids, giving some indication that even that poorly-constrained parameter was indeed measured correctly.

Support for this work was provided by NASA through grant number AR-09196.01 from the Space Telescope Science Institute. All of the data presented in this paper were obtained from the Multimission Archive at the Space Telescope Science Institute (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.