text:
# Stochastic Resonance in Two DimensionalLandau Ginzburg Equation

## Abstract

We study the mechanism of stochastic resonance in a two dimensional Landau Ginzburg equation perturbed by a white noise. We shortly review how to renormalize the equation in order to avoid ultraviolet divergences. Next we show that the renormalization amplifies the effect of the small periodic perturbation in the system. We finally argue that stochastic resonance can be used to highlight the effect of renormalization in spatially extended system with a bistable equilibria.

About twenty years ago, in studying the possible causes of the ice ages sutera1 sutera2 , we introduced the concept of stochastic resonance BSV . This process (hereafter,SR) requires that a physical system has its fundamental symmetry broken in such a away that its phase space is not homogeneous. Thereby, the time spent in different regions of the phase space is not uniformly distributed. If the system is stochastically forced, the probability of large fluctuation can be periodically locked to a small external periodic forcing. The relevant implications of the SR are the generality of its set up and the limited band, in the parameter space, wherein the mechanism occurs. The latter allows for its exploiting in enhancing the signal to noise ratio in a noisy system, while the former makes the concept applicables to a large class of systems. Furthermore, we notice that the mechanism only requires that the system has large deviations properties and, therefore, the restriction to systems stochastically perturbed, is just a case amenable to explicit evaluation of the resonace range in the parameter space. In practice, a stochastic forcing is required only for systems of gradient type (i.e., when a Lyapunov function exists surely) or whenever single point attractors exist.

At the time of the publication above mentioned, we had doubts whether the mathematical process leading to stochastic resonance were a stochastic process at all. Fortunately in f2000  has proved that a wide class of stochastic process obeying to a SR mechanism exists. Therefore we feel confident that, at least in a finite dimensional space, the SR is a well posed mathematical problem. When we consider stochastically perturbed infinite dimensional systems the problem may lead to a not consistent mathematics since the equations for the backward and the foreward probability density function may not exist. Moreover, when a white noise in space and time is considered, the resulting equations are plagued by divergent correlation functions. In a series of papers jona1 ,jona2 ,jona3 , it was proved that these problems could be overcame by the method of the renormalisation theory. On the ground of these rigorous, mathematical results, in the present paper we wish to study the consequences of applying an external periodic (in time) forcing to a simple infinite dimensional stochastic equation.

In particular, we discuss the SR phenomenon in two dimensional Landau Ginzburg equation stochastically perturbed by a white noise:

$$_t\varphi =\nu \mathrm{\Delta }\varphi \frac{V(\varphi )}{\varphi }+ϵ\frac{dW}{dt}$$
(1)
where

$$V(\varphi )=\frac{1}{2}m\varphi ^2\frac{1}{4}g\varphi ^4$$
(2)
and $`\frac{dW}{dt}`$ is a white noise delta correletad both in space and in time. The field $`\varphi `$ is assumed to be real and defined on the torus $`[0,L]X[0,L]`$,$`L=2\pi `$, with periodic boundary conditions. Notice that the previous equation describes mean field theories nearby critical points, therefore, the use of white noise in space and time is mandatory when we do not want to introduce unknown physical scales. Moreover, unless $`m<0`$, the theory has no broken symmetries.

As we mentioned, the phenomenon of SR has been widely studied in literature since it was introduced BPSV ; NIC82 (for a review see rev ). The study of SR in one space dimension, with a broken symmetry, has been presented in BSVJPA . In this case,in fact, the stochastic equation does not show any particular novelty. As we shall see, the extension to more spatial dimension leads to interesting consequences. In particular we will consider their effects on the SR.

As it has been discussed in BJLS  that equation (1) has no meaning in more than one space dimension unless it is renormalized. In D=2, the renomarlization amounts to change any non linear term in (1) by its Wick product, a rather well known procedure in quantum field theory. By using the Wick product, equations (1) and (2) becomes:

$$_t\varphi _\mathrm{\Lambda }=\nu \mathrm{\Delta }\varphi _\mathrm{\Lambda }m\varphi _\mathrm{\Lambda }g\varphi _\mathrm{\Lambda }^3+3gE_m(\mathrm{\Lambda })\varphi _\mathrm{\Lambda }+ϵ\frac{dW_\mathrm{\Lambda }}{dt}$$
(3)
In eq. (3) we have introduced the cutoff field $`\varphi _\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is an ultraviolet cutoff in the Fourier space. Namely by denoting $`\mathrm{\Phi }(\stackrel{}{k})`$ the Fourier transform of $`\varphi (x)`$ we have:

$$\varphi _\mathrm{\Lambda }=_{|\stackrel{}{k}|<\mathrm{\Lambda }}\mathrm{\Phi }(\stackrel{}{k})e^{i\stackrel{}{k}\stackrel{}{x}}𝑑\stackrel{}{k}$$
(4)
and an analogous definition has been used for $`\frac{dW}{dt}`$. In (3) the quantity $`E_m(\mathrm{\Lambda })`$ is the expectation value of second order moment of the field $`Z_\mathrm{\Lambda }`$ satisfying the linear stochastic differential equation:

$$Z_\mathrm{\Lambda }=\nu \mathrm{\Delta }Z_\mathrm{\Lambda }mZ_\mathrm{\Lambda }+ϵ\frac{dW_\mathrm{\Lambda }}{dt}$$
(5)
In two space dimensions one has:

$$E_m(\mathrm{\Lambda })=C\frac{ϵ^2}{\nu }ln[\frac{\nu \mathrm{\Lambda }^2+m}{m}]$$
(6)
where $`C`$ is a constant order one. Finally, let us recall that the quantity $`\varphi _\mathrm{\Lambda }^33E_m(\mathrm{\Lambda })\varphi _\mathrm{\Lambda }`$ is the Wick product of $`\varphi _\mathrm{\Lambda }^3`$.

One can show BJLS  that eq. (3) has a well defined limit for $`\mathrm{\Lambda }`$ going to infinity. Let us remark that this limit is achieved by keeping constant the finite volume. In this case one can show that equation (1) becomes:

$$_t\varphi _\mathrm{\Lambda }=\nu \mathrm{\Delta }\varphi _\mathrm{\Lambda }\frac{V_{eff}(\varphi _\mathrm{\Lambda })}{\varphi _\mathrm{\Lambda }}+ϵ\frac{dW_\mathrm{\Lambda }}{dt}$$
(7)
where:

$$V_{eff}(\varphi )=\frac{1}{2}\rho \varphi ^2\frac{1}{4}g\varphi ^4$$
(8)
The value of $`\rho `$ is determined by the following set of equations:

$$\rho =m+3g[E_m(\mathrm{\Lambda })E_\alpha (\mathrm{\Lambda })]$$
(9)
$$\alpha =\frac{^2V_{eff}}{\varphi ^2}at\varphi =\varphi _e$$
(10)
where $`\varphi _e`$ is the equilibria of the effective potential, i.e. $`\varphi _e`$ is determined by solving the equation:

$$\frac{V_{eff}}{\varphi }=0$$
(11)
The set of equations (9,10,11) have a clear physical meaning BJLS , shortly reviewed as follows. First of all, let us note that $`\rho `$ is independent of $`\mathrm{\Lambda }`$ in the limit $`\mathrm{\Lambda }\mathrm{}`$ because of eq. (6). Eq. (9) is a redefinition of the Wick product in terms of a new constant $`\alpha `$, which is defined in equation (10) as the curvature of the effective potential at one of its equilibrium values, namely $`\varphi _e`$. The physical idea of eq. (10) is that, the counter term introduced in order to renormalize the Landau Ginzburg equation, produces a double well potential. The position of the new minimum, determined by $`\rho `$, must be consistent with the fluctuations around the minimum itself.

Using equations (9,10,11) one gets a nonlinear equation for the variable $`x=\frac{\alpha }{m}`$, namely

$$x=2+2Klogx$$
(12)
where $`K=\frac{3gCϵ^2}{m\nu }`$. Let us remark, for future purpose, that it is possible to show that a solution of eq. (12) always satisfies $`x>2K`$. As discussed in BJLS , the effective potential describes rather well the statistical properties of the space averge field of eq.(1).

In general we can consider the variables $`\mathrm{\Phi }(l)`$ defined as the average in space of $`\varphi (x)`$ on a box of side $`l`$. For $`l<L=2\pi `$, the statistical properties of $`\mathrm{\Phi }(l)`$ are no longer described by the same value of $`\rho `$ because of the increasing effect, at small scales, of the stochastic fluctuations. In particular, for $`l0`$, one finds that the statistical properties of $`\mathrm{\Phi }(l)`$ are described by a single well potential centered around $`\mathrm{\Phi }(l)=0`$.

We want now to discuss the solutions of eq. (1) when a small periodic perturbation is added to the system and in particular we want to discuss the effect of renormalization of the mechanism on the stochastic resonance.

To this aim, following the original discussion given in BSV  we consider the case of a time independent constant $`A`$ added to the r.h.s. of (1). The crucial point is the computation of the effective potential $`V_{eff}`$ previously introduced. Equation (9) are unchanged while equations (10,11) become:

$$\alpha =\rho +3g\varphi _{e}^{}{}_{}{}^{2}$$
(13)
$$\rho \varphi _eg\varphi _{e}^{}{}_{}{}^{3}+A=0$$
(14)
The effect of the constant perturbation is to change the equilibrium solution of the effective potential and therefore to change the fluctuations around the equilibrium solution, which in turn change the value of $`\rho `$. The latter effect is due only to renormalization.

In order to understand qualitatively the contribution due to renormalimation, one can compute a perturbative solution of the set of equations (9,13,14) in power of $`A`$. At first order one obtains:

$$\rho _{1,2}=\rho _0+\frac{6g\varphi _{e0}H}{12H}\frac{A}{2\rho _0}$$
(15)
$$\varphi _{e1,e2}=\varphi _{e0}+\frac{1+H}{12H}\frac{A}{2\rho _0}$$
(16)
where $`H=\frac{3gCϵ^2}{\nu \alpha _0}`$ ( $`\alpha _0=2\rho _0`$) and all quantities with the index $`0`$ refer to values computed for $`A=0`$. The previous definition of $`H`$ implies $`H=\frac{K}{x}<\frac{1}{2}`$ because of the discussion on the solution of eq. (12). Let us remark that we obtain two different values of $`\rho `$ and $`\varphi _e`$ because two different values of $`\varphi _{e0}`$ are possible, one negative and the other one positive. Thus there exists two possible effective potential, each one describing the two well respectively.

By using equations (15) and (16) is now possible to compute the effective potential difference between the two stable solutions and the unstable solution. One gets:

$`\mathrm{\Delta }V_{1,2}={\displaystyle \frac{\rho ^2}{4g}}+A_R\varphi _{e0}`$ (17)
$`A_R=A(1+{\displaystyle \frac{3H}{2(12H)}})`$ (18)
In the limit of $`ϵ0`$ we have $`H0`$ and the effect of renormalization disappears. For finite value of the noise and for $`g`$ large enough the renormalization amplifies the asymmetry of the double well potential. Therefore, we should expect that, in two dimensions Landau Ginzburg equation, the mechanism of SR acts for smaller value of the forcing amplitude with respect to the one estimated by looking at the equilibrium probability distribution with no periodic forcing.

In order to test the renormalization effect on SR, we have performed a numerical simulation on square lattice of 16x16 points. The parameters of our numerical simulation are $`g=220`$, $`\nu =0.1`$, $`m=0.1`$, $`ϵ=0.1`$. On the lattice, the quantity $`E_m(\mathrm{\Lambda })`$, defined in eq. (6), is given by:

$$E_m(a)=\mathrm{\Sigma }_{p,q}\frac{a^2ϵ^2}{2S(m,\nu ,q,p)}$$
(19)
where

$$S(m,\nu ,q,p)=ma^2+4\nu 2\nu cos(2\pi q/N)2\nu cos(2\pi p/N)$$

In equation (19 we introduce the lattice mesh $`a=2\pi /N`$. It follows that $`\mathrm{\Lambda }=a^1`$.

For different values of $`A`$, by applying Newton method, we have computed the solution of equations (9,13,14), using (19). It turns out that, with our numerical parameters, $`\rho _0=4.9.`$ and $`A_R2A`$, i.e, the external forcing is amplified by almost a factor 2.

We have performed a long time integration of (1) with $`A=0`$. The equilibrium probability distribution of $`\mathrm{\Phi }(L)`$ has been found to be bimodal with maxima centered at $`\pm \sqrt{\frac{\rho _0}{g}}`$ and the average transition time $`\tau `$ between the two equilibrium $`\tau =9`$ in time units of $`m^1`$. Next, using the analytical theory of SR discussed in BSV , we have computed, by taking into account the amplification of the renormalization, the value of $`A`$ for which the system should show a clear SR. We have found $`A=0.1`$. Finally we have applied a periodic perturbation on the r.h.s. of eq. (1) with amplitude $`A=0.1`$ and period equal to $`2\tau =18`$, i.e. two times the average transition time. We have performed a long time integration corresponding to 2048 time units. In figure (1) we show the power spectrum of the average field $`\mathrm{\Phi }(L)`$. As one can clearly seen, there is an extremely well defined peak at frequency 113 which corresponds to a period of 18 time units.

In order to show the amplification effect on SR due to renormalization, we have numerically integrated the zero dimensional stochastic differential equation:

$$\frac{dx}{dt}=\rho _0xgx^3+Acos(\frac{2t\pi }{18})+\sigma \frac{dW}{dt}$$
(20)
where $`A=0.1`$ as before and the variance $`\sigma `$ of the white noise has been tuned in order to reproduce an average transition time of $`9`$ ($`\sigma =0.0165`$) for $`A=0`$. With this choice of $`\sigma `$ we have been able to reproduce the equilibrium probability distribuiont of $`\mathrm{\Phi }(L)`$ at $`A=0`$ for the 2D Landau Ginzburg equation. As for the previous case, we have integrated eq. (20) for $`2048`$ time units. In figure (2) we show the power spectrum of $`x`$. At variance with figure (1), only a small effect of the periodic forcing is felt by the system, i.e. the system does not exhibits any SR.

The comparison between figures (1) and (2), indicates clearly that the renormalization effect is acting as an amplifier of the external forcing, as previously discussed in the framework of a perturbation theory.

It is quite interesting to compute the statistics of the transition times between the minima at different scales as described by the variables $`\mathrm{\Phi }(l)`$. In absence of forcing the probability distribution of the transition time for $`\mathrm{\Phi }(L)`$ is exponential BJLS . In the stochastic resonance the probability distribution of the transition time for $`\mathrm{\Phi }(L)`$ should show a well defined maximum at 9 time units. On the other hand, if we consider the transition time for $`\mathrm{\Phi }(l)`$ with $`lL`$, then we expect that the forcing is not able to produce any SR. This is due to the fact that the parameter describing the effective potential for small $`l`$ are changed and, in particular, the average transition time between the minima becomes much smaller BJLS . This effect is indeed observed in our numerical simulation. In fig. (2) we show the probability distribution of the transition time for $`\mathrm{\Phi }(L)`$ and $`\mathrm{\Phi }(\frac{L}{16})`$. As predicted, for the small scale, the probability distribution of the transition time is exponential.

We can therefore reach the conclusion that, becuase of the renormalization, the SR mechanism becomes scale depedent, a feature not previously observed.

In this letter we have shown that renormalization can amplify the SR in two dimensional Landau Ginzburg equation. We want to remark that this effect can also be used to test, in a given experimental situation, whether or not a renormalization mechanism is acting in the system. In this case SR can be used as a tool to measure quantitatively the effect, if any, of renormalization in spatially exnteded systems. We remark that the present effect may be not restricted to system of gradient type and, in principle, could be detected in spatially extended deterministic systems.

Acknowledgements We thank Prof. Jona Lasinio for useful discussions.


format: text

text:
# Correlations between the WMAP and MAXIMA Cosmic Microwave Background anisotropy maps

## 1. Introduction

Temperature fluctuations in the cosmic microwave background (CMB) encode a vast amount of cosmological information about our universe. CMB photons released from the primordial plasma at the time of recombination approximately 380,000 years after the Big Bang, provide thus a picture of the universe in its infancy only somewhat modified by low-redshift effects such as reionization. Recently WMAP produced a 13 full sky measurement of the CMB temperature anisotropy (Bennett et al., 2003a). This map has been used in conjunction with other CMB and cosmological data to constrain a number of cosmological parameters to unprecedented accuracy (Spergel et al., 2003). Previous to WMAP a number of experiments produced high quality maps of CMB temperature anisotropy. These included both balloon borne bolometric experiments such as BOOMERANG (de Bernardis et al., 2000; ruhl/etal:2002, Ruhl et al.(2002), MAXIMA-I (Hanany et al., 2000; Lee et al., 2001), and ARCHEOPS (Benoit et al., 2003), and ground based interferometric experiments, CBI (Padin et al., 2001; Mason et al., 2003), DASI (Halverson et al., 2002), and VSA (Grainge et al., 2002). Tight constraints were placed on cosmological parameters from these experiments as well (e.g. Jaffe et al. (2001); Netterfield et al. (2002); Stompor et al. (2001); Abroe et al. (2002); Pryke et al. (2002)).

Given the longer observations and higher sensitivities of recent CMB experiments the theoretical analysis is more likely to be limited by systematic rather than statistical errors. It is therefore important to ensure that systematic errors in these experiments are sub-dominant compared to statistical errors. A comparison of the power spectra from the experiments can provide some confidence that systematic errors are not dominant. For instance the power spectra of WMAP, MAXIMA-I, and MAXIMA-II are shown in Figure 1. Note that no calibration adjustments have been made to the data. For experiments which observe overlapping parts of the sky the close agreement of the power spectra does not necessarily imply that the spatial fluctuations detected by the experiments are identical. In such cases, and particularly when the experiments have similar angular resolution, a more stringent test for systematic errors is to cross correlate the temperature fluctuations of one map with the fluctuations in the other. Positive correlations between temperature anisotropy maps would also enhance the confidence in the reconstruction of the spatial pattern of the CMB. Some difficulties arise if the two experiments under consideration have different pixel resolutions and beam profiles, which is usually the case. In that case a straightforward pixel to pixel comparison is no longer accurate because the CMB signal contained in corresponding pixels is not the same, and a more elaborate technique needs to be employed.

In this paper we use the cross-spectrum as a technique to cross-correlate the maps of WMAP and MAXIMA. This is the first reported cross-correlation of a CMB map with the map of WMAP, and the first release of data from the MAXIMA-II flight.

Other comparisons of CMB anisotropy maps from different experiments have been performed in the past (e.g. Ganga et al. (1993), Tenerife and DMR (Linweaver et al., 1995), MSAM and Saskatoon (Knox et al., 1998), and QMAP and Saskatoon (Xu, Tegmark, & de Oliveira-Costa, 2002)). However, none of these analyses used the cross spectrum as a technique for quantifying the amount correlation between data sets.

This paper is organized as follows: in Section 2 we give details of the cross spectrum and its use for quantifying the amount of correlations between two CMB anisotropy maps at various angular scales. In Section 3 we discuss the maps used in our analysis. Results of the calculation of the power and cross-spectra for the maps, and an analysis of the difference maps are given in Section 4. We discuss the results in Section 5, and summarize in Section 6.

## 2. Method

Perhaps the simplest way to compare quantitatively two maps is to calculate a $`\chi ^2`$ for their difference. This statistic has the limiting feature that only overlapping sections of the maps can be included, and the maps must have identical pixelizations. Additionally, if the maps have different beam window functions, then the signal in corresponding pixels is different and the distribution of the $`\chi ^2`$ statistic would no longer be $`\chi ^2`$ distributed. The statistical interpretation of the $`\chi ^2`$ value would therefore be difficult.

There are a number of statistics which can be calculated between two CMB maps to quantify their consistency, e.g. the linear and rank correlation coefficients (Press et al., 1992). Unfortunately, these correlation coefficients also suffer from the same difficulties as the $`\chi ^2`$ statistic.

Several authors derived statistics which do take into account partial overlapping maps, different pixelizations and beam profiles. Knox et al. (1998) derive both Bayesian and frequentist techniques for correlating CMB maps, and apply these statistics to data from the MSAM92, MSAM94, and Saskatoon experiments. They advocate the calculation of the “contamination parameter”, which gives the probability distribution for the magnitude of a signal that is not common between the two data sets under consideration. A low value for the contamination parameter implies that the data sets are consistent. Tegmark (1999) defines a “null-buster” statistic which gives the number of $`\sigma `$ between the difference of two maps and a hypothesis of pure noise. These statistics can be used both as an internal consistency check between detectors for the same experiment (Stompor et al., 2003) and to compare maps from different experiments.

In this paper we use the cross-spectrum to quantify the level of correlations between the data of WMAP and MAXIMA . The cross-spectrum is the spherical harmonic transform of the real space correlation function for the two maps. Whereas other statistics condense the information about correlations into a single number, and therefore result in some loss of information, the cross-spectrum retains more information by analyzing the correlations as a function of angular scale.

We now discuss the method for estimating the cross-spectrum from CMB temperature anisotropy maps. Consider two maps of the CMB called $`T^{(\mathrm{𝟏})}`$ and $`T^{(\mathrm{𝟐})}`$, respectively. Then

$`\mathrm{\Delta }T_i^{(\mathrm{𝟏})}`$ $`=`$ $`s_i^{(\mathrm{𝟏})}+n_i^{(\mathrm{𝟏})}`$ (1)
$`\mathrm{\Delta }T_i^{(\mathrm{𝟐})}`$ $`=`$ $`s_i^{(\mathrm{𝟐})}+n_i^{(\mathrm{𝟐})},`$ (2)
where $`i`$ is a pixel index, $`s_i^{(\mathrm{𝟏})}`$ and $`s_i^{(\mathrm{𝟐})}`$ are the CMB signal in pixel $`i`$, and $`n_i^{(\mathrm{𝟏})}`$ and $`n_i^{(\mathrm{𝟐})}`$ are the pixel noise for the first and second map, respectively. Let the number of pixels in the first and second maps be $`N_\mathrm{p}^{(\mathrm{𝟏})}`$ and $`N_\mathrm{p}^{(\mathrm{𝟐})}`$, respectively. We write the data vector in pixel space as

$`d`$ $`=`$ $`\left[\begin{array}{c}\mathrm{\Delta }T^{(\mathrm{𝟏})}\\ \mathrm{\Delta }T^{(\mathrm{𝟐})}\end{array}\right],`$ (5)
where $`d`$ is now a column vector of length $`N_\mathrm{p}^{(\mathrm{𝟏})}+N_\mathrm{p}^{(\mathrm{𝟐})}`$. Assuming the signal and noise within each experiment are uncorrelated and that the noise between experiments is uncorrelated, and using Equations 1 and 2 we find that

$`dd^T`$ $`=`$ $`M\left[\begin{array}{cc}S^{(\mathrm{𝟏})}+𝒩^{(\mathrm{𝟏})}& S^{(𝐂)}\\ S^{(𝐂)T}& S^{(\mathrm{𝟐})}+𝒩^{(\mathrm{𝟐})}\end{array}\right],`$ (8)
where $`𝒩^{(\mathrm{𝟏})}=n^{(\mathrm{𝟏})}n^{(\mathrm{𝟏})T}`$ and $`𝒩^{(\mathrm{𝟐})}=n^{(\mathrm{𝟐})}n^{(\mathrm{𝟐})T}`$ are the pixel noise covariance matrices for the first and second experiment, respectively. The quantities $`S^{(\mathrm{𝟏})}`$, $`S^{(\mathrm{𝟐})}`$, and $`S^{(𝐂)}`$ are the CMB signal covariance matrices, which can be written as

$`S_{ij}^{(\mathrm{𝟏})}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{}}{\displaystyle \frac{2\mathrm{}+1}{4\pi }}C_{\mathrm{}}^{(\mathrm{𝟏})}B_{\mathrm{}}^{(\mathrm{𝟏})2}P_{\mathrm{}}(\mathrm{cos}\theta _{ij})`$ (9)
$`S_{ij}^{(\mathrm{𝟐})}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{}}{\displaystyle \frac{2\mathrm{}+1}{4\pi }}C_{\mathrm{}}^{(\mathrm{𝟐})}B_{\mathrm{}}^{(\mathrm{𝟐})2}P_{\mathrm{}}(\mathrm{cos}\theta _{ij})`$ (10)
$`S_{ij}^{(𝐂)}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{}}{\displaystyle \frac{2\mathrm{}+1}{4\pi }}C_{\mathrm{}}^{(𝐂)}B_{\mathrm{}}^{(\mathrm{𝟏})}B_{\mathrm{}}^{(\mathrm{𝟐})}P_{\mathrm{}}(\mathrm{cos}\theta _{ij}),`$ (11)
where $`B_{\mathrm{}}^{(\mathrm{𝟏})}`$ and $`B_{\mathrm{}}^{(\mathrm{𝟐})}`$ are the beam profiles for the first and second experiment, $`P_{\mathrm{}}`$ are Legendre polynomials, and $`\theta _{ij}`$ is the angle between pixels $`i`$ and $`j`$. The auto-spectra $`C_{\mathrm{}}^{(\mathrm{𝟏})}`$ and $`C_{\mathrm{}}^{(\mathrm{𝟐})}`$ are commonly called the power spectra for the first and second map, respectively, and $`C_{\mathrm{}}^{(𝐂)}`$ is defined as the cross-spectrum. In the zero noise case the cross-spectrum is limited by the requirement that $`C_{\mathrm{}}^{(𝐂)}\sqrt{C_{\mathrm{}}^{(\mathrm{𝟏})}C_{\mathrm{}}^{(\mathrm{𝟐})}}`$. Otherwise $`M`$ would have negative eigenvalues and therefore be unphysical. Up to an irrelevant additive constant, the likelihood for $`d`$ is

$`2\mathrm{ln}`$ $`=`$ $`(N_p^{(\mathrm{𝟏})}+N_p^{(\mathrm{𝟐})})\mathrm{ln}|M|+d^TM^1d.`$ (12)
We maximize the likelihood $``$ as a function of the two auto-spectra and the cross-spectrum. If two CMB maps have a high degree of correlation, then we expect the cross-spectrum to resemble the auto-spectra of both maps. If the signal in the two maps is not correlated the cross-spectrum should be consistent with zero at all angular scales.

To maximize $``$ we adopt a Newton-Raphson technique for finding the zero of the first derivative of the log likelihood function. We use a technique similar to the one described by Bond, Jaffe, & Knox (1998) except that we use the full curvature matrix instead of the Fisher matrix when calculating the steps for convergence of the Newton-Raphson algorithm (Hobson & Maisinger, 2002). We found that use of the Fisher matrix gave non-positive definite pixel-correlation, presumably because the likelihood function of the three spectra together has a complicated structure.

The cross-spectrum method involves estimating all three power spectra simultaneously and gives rise to correlations between the different spectra. Therefore an auto-spectrum estimated for any one experiment alone may differ from the auto-spectrum calculated when estimating the cross spectrum. For the case of a full sky coverage one can show that the level of correlations between the spectra is proportional to the amplitude of the cross-spectrum, that is to the amount of common sky signal (Kamionkowski et al., 1996). We therefore expect a high level of correlations between the spectra if they share the same sky signal at low $`\mathrm{}`$’s where the contribution of instrument noise is smaller compared to sample variance.

The cross-spectrum is a powerful technique for comparing two CMB anisotropy maps because it accounts for different beam shapes, pixel resolutions and sky coverage in a simple and straightforward way. For example, one can compute the cross-spectrum for two maps that do not overlap at all. In such a case the results would not be sensitive to correlations on small angular scales. Also, though our formalism describes estimating the cross-spectrum for two maps, a further generalization could be made to estimating correlations between three or more maps simultaneously. However, given numerical subtleties which are subsequently discussed in Section 5.2, we found it prudent, while equally convincing, to estimate the correlations for only pairs of maps.

The power spectra estimated from the WMAP maps are actually the cross-spectra from different detectors from the same frequency band (Hinshaw et al., 2003). Our approach is distinct from the one used by the WMAP team in that they use a frequentist approach (Hivon et al., 2002) to estimate the cross-spectra between various detectors of the same frequency band independently from the auto-spectra. The approach we use in this paper is entirely Bayesian.

## 3. The Maps

We cross-correlate CMB temperature anisotropy maps from the WMAP, MAXIMA-I and MAXIMA-II experiments. The WMAP map used is the W-band (93 GHz) foreground-cleaned map, and only the portion that overlaps with the MAXIMA-I field. By analyzing
WMAP ’s derived foreground maps we find that this portion of the W-band map is free of point sources, and contains a negligible amount of dust, free-free, and synchrotron emission111http://lambda.gsfc.nasa.gov/product/map/m\_products.cfm. The WMAP data at 93 GHz has similar angular resolution and is closest in frequency to the MAXIMA data. This map is pixelized using the HEALPix222http://www.eso.org/science/healpix/ (Gorski, Hivon, & Wandelt, 1999) pixelization in Celestial coordinates with nside$`=512`$, and contains 7,926 pixels. An nside=512 corresponds roughly to a pixel with a width of 7. Following Bennett et al. (2003a) we assume that there are no noise correlations between pixels and that the beam pattern is non-Gaussian, azimuthally symmetric and has FWHM of 13.2. We computed the beam window function by taking the weighted average of the individual beams for each of the four detectors used to form the final W-band map. The resulting beam profile in $`\mathrm{}`$ space is shown in Figure 2. The Wiener filtered version of this map is shown in the left panel in Figure 3; we use the raw data in the analysis. The Wiener filtering is performed with the corresponding best fit models for each map. A more detailed discussion of the WMAP maps is given in Bennett et al. (2003a) and of algorithms used in their computations in Hinshaw et al. (2003).

The MAXIMA map-making procedure is described exhaustively in Stompor et al. (2002). The MAXIMA-I map we use in this analysis is the 8 arcminute version of the data published by Hanany et al. (2000), which covers a larger area of the sky and has a coarser resolution compared to the data published by Lee et al. (2001). There are a total of $`5,972`$ pixels, and this map covers $`100`$ $`\mathrm{deg}^2`$ on the sky. It is a combination of three $`150`$ GHz photometers and one 240 GHz photometer. The beams for each photometer have a FWHM of $`10^{}`$ (Hanany et al., 2000) and the effective window function is calculated using the technique described in Wu et al. (2001). The MAXIMA-I map shown in Figure 3 is Wiener filtered and smoothed to a WMAP resolution. The raw version of the map is used for all quantitative analyses. Figure 3 also shows the difference of the MAXIMA-I and WMAP Wiener filtered maps. The pattern of temperature fluctuations, which is similar in both maps, disappears in the difference map.

The MAXIMA-II map comes from a flight of the MAXIMA payload that took place on 1999 June 17. We use the data from four photometers at 150 GHz, and only the portion of the MAXIMA-II map that overlaps the map of MAXIMA-I. More details about the MAXIMA-II flight, data and maps are given in Rabii et al. (2003, in preparation), and Stompor et al. (2003). The MAXIMA-II map is pixelized using an 8 square pixelization in celestial coordinates, contains 2,757 pixels, and covers $``$ 50 $`\mathrm{deg}^2`$ on the sky. The beam profile for this map is $`10^{}`$ FWHM, and again computed using the techniques described in Wu et al. (2001). The MAXIMA-II power spectrum shown in Figure 1 has 10 bins of $`\mathrm{\Delta }\mathrm{}=75`$, extending over the $`\mathrm{}`$ range $`35\mathrm{}785`$. Figure 4 shows the overlap region of the MAXIMA-I and MAXIMA-II maps and the difference map. Identical temperature fluctuations that are apparent in each of the maps disappear in the difference map.

## 4. Results

The auto- and cross-spectra for all combinations of the WMAP and MAXIMA maps are shown in Figure 5. The error bars on the spectra are the square root of the curvature of the likelihood function about the maximum likelihood parameter value. In all cases we compute the spectra in bins of width $`\mathrm{\Delta }\mathrm{}=75`$, over the interval $`111\mathrm{}710`$, and marginalize over all modes $`\mathrm{}110`$ and $`\mathrm{}711`$. The appropriate pixel window functions for each map were convolved with the beam functions in the analysis. We found that the cross-spectrum estimator did not converge when the initial bin was split in two, and this is further discussed in Section 5. In all cases the cross-spectra are consistent with the auto-spectra giving strong evidence for a correlation between the maps.

We also compute the power spectrum of the difference maps for all three pairs of maps using bins of $`\mathrm{\Delta }\mathrm{}=75`$ over the range $`35\mathrm{}785`$. The WMAP window function is used when computing the MAXIMA-I -WMAP and MAXIMA-II -WMAP difference spectra. The expected residual power resulting from different beam profiles is maximum at the bin centered at $`\mathrm{}300`$, and is approximately equal to the 1$`\sigma `$ error bar of the MAXIMA-I /WMAP difference power spectrum. The effect is less than 1$`\sigma `$ for all remaining bins. We use the MAXIMA-I window function when computing the MAXIMA-I -MAXIMA-II difference spectrum. The results are shown in Figure 6. Of the 30 band power estimates for the difference maps, 28 are within 1$`\sigma `$ of zero power.

To further quantify the level of correlation between the maps we use a $`\chi ^2`$ statistic to reject the hypothesis that the maps are uncorrelated. We write our statistic as

$`\chi ^2={\displaystyle \underset{BB^{}}{}}C_B^{(𝐂)}F_{BB^{}}C_B^{}^{(𝐂)}`$ (13)
where the sum is over band power estimates, and $`F`$ is the Fisher matrix for the cross-spectrum. Because the auto-spectra and cross-spectrum are estimated simultaneously, we marginalize over the auto spectra when calculating the $`\chi ^2`$.

To test the null hypothesis we choose a statistical significance $`\alpha =1\times 10^8`$. If $`\chi ^2`$ is greater than the critical value $`53`$, then the probability that the null hypothesis is true is less than $`1\times 10^8`$. The results are summarized in Table 1. In all cases $`\chi ^2`$ is significantly larger than the critical value giving an essential certainty that the no-correlation hypothesis is false. Note that assuming the $`\chi ^2`$ of Equation 13 is $`\chi ^2`$ distributed is equivalent to assuming that the $`C_B^{(𝐂)}`$ are Gaussian distributed, which is an approximation.

We also compute the $`\chi ^2`$ of the difference spectra shown in Figure 6 with the null spectra to determine how consistent these are with no fluctuations in the difference maps. The results are shown in Table 2. The 10 power spectrum bins computed from the MAXIMA-I /WMAP and MAXIMA-I /MAXIMA-II difference maps have a $`\chi ^2`$ of 7.5 and 8.3, respectively, with the null spectrum. The MAXIMA-II /WMAP difference map gives a $`\chi ^2`$ of 17.2. There is a 7% chance of getting $`\chi ^217.2`$ for 10 DoF. Overall there is a good fit to the null spectrum model, which implies that differencing the overlap section of the maps removes the sky signal and is consistent with noise.

## 5. Discussion

### 5.1. Auto- and Cross-Spectra

The auto- and cross-spectra of the different data sets agree with each other to within $`1\sigma `$ over almost all $`\mathrm{}`$ bins giving evidence that at each angular scale all experiments are detecting the same spatial fluctuations on the sky. All auto- and cross-spectra show the first acoustic peak in the power spectrum and then a level of power that is consistent with subsequent peaks. These results are consistent with standard inflationary $`\mathrm{\Lambda }`$CDM models.

Auto-spectra of the overlap section of the WMAP data give increased error bars at $`\mathrm{}486`$ because of the limited sky coverage of the overlap regions, and because of the beam profile of the W-band map. We find that the beam pattern alone causes the WMAP auto-spectrum error bars in the bins $`\mathrm{}`$=$`\{\{486,560\},\{561,635\},\{636,710\}\}`$ to be 2-3 larger than those for MAXIMA-I or MAXIMA-II. Negative power was found in the bin $`\mathrm{}`$=$`\{486,560\}`$ for the WMAP auto-spectra (see the top and bottom panel of 5). There is no requirement that the auto-spectrum be positive in our estimation method.

A comparison of the auto-spectra shown in Figure 5 reveals that there is a difference between band power estimates for the same dataset. This difference arises because the computation of cross-spectra involves estimating both auto- and cross-spectra simultaneously, giving rise to correlations between the different spectra. The fractional changes in power averaged over bins are 3%, 6%, and 15% for the WMAP, MAXIMA-I, and MAXIMA-II data sets, respectively. If the likelihood distribution of the band powers were strictly Gaussian, then the maximum likelihood estimates would be the same regardless of the correlations between spectra. However, the likelihood as a function of auto- and cross-spectra is somewhat non-Gaussian (see Equation 12, Bond, Jaffe, & Knox (1998)), so the correlations do effect the band power estimates. The fact that the changes between estimates are small suggests, however, that the distributions are close to Gaussian.

We note that the bin-powers are correlated at the 77% level or higher at the lowest $`\mathrm{}`$ bin. The correlation decreases to less than 10% at the highest $`\mathrm{}`$ bin for correlations between the auto-spectra and to 20% - 50% for correlations between the auto- and cross-spectra. These results are in broad agreement with expectations given the high amplitude of the cross-spectrum.

### 5.2. Computational Issues

The strong correlation between the different data sets leads to some computational difficulties when attempting to find the maximum likelihood auto- and cross-power spectra. As discussed in Section 2 the cross-spectrum is limited by the requirement that $`C_{\mathrm{}}^{(𝐂)}\sqrt{C_{\mathrm{}}^{(\mathrm{𝟏})}C_{\mathrm{}}^{(\mathrm{𝟐})}}`$. We find that using a quadratic estimator (Fisher matrix) method for calculating the Newton-Raphson step leads to a mis-estimate of the step for $`\delta C_{\mathrm{}}^{(𝐂)}`$ when starting with a guess significantly far away from the peak in likelihood space. This consequently results in a step which leads to a non-positive definite pixel covariance matrix. For example, an initial guess of a null spectrum leads to an unphysical pixel covariance matrix in all three cases we are considering. This is remedied by using the curvature matrix to compute $`\delta C_{\mathrm{}}^{(𝐂)}`$, as discussed in Section 2. Once the parameter values become sufficiently close to the maximum likelihood values either the curvature matrix or fisher matrix can be used to find the maximum likelihood parameters. Both techniques converge to the same set of parameters for all three analyses.

The power spectra shown in Figure 5 have been calculated with a broad initial bin at $`\mathrm{}=\{2,110\}`$. This is because the Newton-Raphson likelihood maximization technique did not converge if this low $`\mathrm{}`$ bin was split to two. Some binning structures would cause negative eigenvalues in the curvature matrix of the parameters or steps in parameter space that would lead to a non-positive definite pixel covariance matrix, both of which are unphysical.

We carried out simulations and found a similar phenomenon. The cross-spectra of uncorrelated maps or maps with a small value for the ratio of expected cross-spectrum to auto-spectrum converged to the expected answer. The calculation of the cross-spectrum also converged with simulated maps that had perfect correlation (i.e. the same map with different noise realizations) and a broad first bin with $`\mathrm{}=\{2,110\}`$. However it did not converge with simulated maps that had perfect correlation and two bins between $`\mathrm{}`$ of $`2`$ and $`110`$. Therefore, we attribute the computation problems encountered as a limitation in the method used for computing the cross-spectrum and not a feature in any of the data sets considered in this analysis.

### 5.3. Foregrounds and Systematic Errors

The cross spectrum and difference spectrum analyses conclusively demonstrate that all three experiments have mapped the same temperature fluctuations on the sky. However, these analyses are not sensitive to whether the shared fluctuations are CMB in origin or the result of foreground contamination or a shared systematic effect.

A careful foreground analysis was carried out by both the WMAP and MAXIMA teams. It was shown that the MAXIMA-I region of the sky at 150-240 GHz contains a negligible amount galactic contamination and that it has no detectable point sources (Hanany et al., 2000; Jaffe et al., 2003).

A detailed analysis of the foreground sources in the WMAP data is presented in Bennett et al. (2003b). Although we use WMAP ’s foreground-cleaned map,
even WMAP ’s foreground maps have negligible amount of contamination in the MAXIMA-I region. The RMS fluctuations in the WMAP 93 GHz dust, synchrotron, and free-free maps are lower than the corresponding WMAP CMB map by a factor of $`17`$, $`42`$, and $`560`$, respectively. Also, the WMAP team find no point sources in the MAXIMA-I region of the sky.

It is unlikely that WMAP and MAXIMA share systematic errors. We therefore conclude that the common signal in the WMAP and MAXIMA data is the cosmic microwave background radiation. Since the cross-spectra agree with the auto-spectra we conclude that within the signal-to-noise ratio of the tests systematic errors in the data are smaller compared to statistical errors.

## 6. Summary

We have presented a Bayesian method for estimating the cross-spectrum between two CMB temperature anisotropy maps. The method is advantageous for correlating maps because it does not require the maps to have perfect overlap, identical beam shapes or pixelizations. Using this formalism we found a high degree of correlation between the maps from MAXIMA-I, WMAP, and MAXIMA-II ; in all cases the null hypothesis is rejected with a probability higher than $`110^8`$. Additionally, we computed the power spectrum of the difference maps for all combinations of the three data sets considered, and found that in each case the spectra were consistent with the null spectrum.

The results show conclusively that the temperature fluctuations detected by each of the MAXIMA-I, WMAP, and MAXIMA-II experiments are reproduced by these experiments, in overlapping regions of the sky. The close agreement of the fluctuations detected by these experiments shows that current CMB experiments are now beginning to provide us with high precision images of the true microwave sky.

All computations for this analysis were performed at the University of Minnesota Supercomputing Institute in Minneapolis, Minnesota, and at the National Research Scientific Computing Center in Berkeley, California, which is supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC03-76SF00098. We acknowledge the use of the HEALPix software package. MEA and RS acknowledge support from NASA grant no. S-92548-F.


text:
# Coherent X-ray Diffractive Imaging; applications and limitations

## I Introduction

Three ideas, developed over half a century, have now converged to provide a working solution to the non-crystallographic phase problem. In this paper we outline some applications and our development of this solution, which provides a method for lensless, diffraction-limited, aberration-free X-ray imaging of nano-objects in three-dimensions at high resolution. We suggest the acronym CXDI (Coherent X-ray Diffractive Imaging) for this field. We also discuss the fundamental limitations on resolution set by radiation damage, and some new approaches to this problem based on femtosecond diffraction.

The ideas start with Sayre’s 1952 observation that Bragg diffraction undersamples diffracted intensity relative to Shannon’s theorem sayre:1952 . Secondly, the development of iterative algorithms with feedback in the early nineteen-eighties produced a remarkably successful optimization method capable of extracting phase information from adequately sampled intensity data fienup:1982 . Finally, the important theoretical insight that these iterations may be viewed as Bregman Projections in Hilbert space has allowed theoreticians to analyze and improve on the basic Fienup algorithm elser:2003 . A parallel development was the real space algorithm for “phase recovery” in crystallography szoke:1997  and non-periodic diffraction hauriege:2003 .

Experimental work started for X-rays with the images of lithographed structures reconstructed from their soft-X-ray diffraction patterns by Miao et al. in 1999 miao:1999 . More recently we have seen higher resolution X-ray imaging miao:2002 , 3D imaging miao:2002 ; williams:2003 , the introduction of algorithms which do not require a conventional lower-resolution image to provide the boundary (support) of the object he:acta:2003 ; he:prb:2003 ; marchesini:2003  and the striking atomic-resolution image of a single nanotube reconstructed from an experimental electron micro-diffraction pattern (an example of the more general CDI method) zuo:2003 .

The rapid growth of nanoscience (described recently as “the next industrial revolution”) has produced an urgent need for techniques capable of revealing the internal structure, in three dimensions, of inorganic nanostructures and large molecules which cannot be crystallized (such as the membrane proteins of vital importance for drug delivery). Scanning probe methods are limited to surface structures, and the electron microscope can provide atomic resolution images of projections of crystalline materials in thicknesses up to about 50nm, or tomography of macromolecular assemblies and inorganics at lower resolution. No technique at present can provide three-dimensional imaging at nanometer resolution of the interior of particles in the micron size range. The development of such a method would have a decisive impact on several fields of science, and would be a timely enabling technology for the Molecular Foundary at Lawrence Berkeley National Laboratory, the Linac Coherent Light Source (LCLS) at Stanford and other nanoscience user facility initiatives. Briefly, we foresee initial applications for CXDI (in the 0.5–4kV X-ray range) in materials science as follows: 1. The visualization of the internal labyrinth structure of the new mesoporous framework structures (eg. glassy foams, now finding uses for molecular sieves and hydrogen storage; 2. Imaging the complex tangles of dislocation lines which are responsible for work-hardening; 3. Imaging the cavities within duplex steels, responsible for their very high uniform extension; 4. 3D imaging defect structures in magnetic multilayers; 5. The tomographic imaging of misfit dislocations at interfaces, free of the thin-film elastic relaxation processes which distort the images obtained by transmission electron microscopy; 6. Imaging of the three-dimensional arrangement of Orowan dislocation loops which, by entanglement with particles, provide the dispersion-hardening of copper alloys; 7. The imaging of precipitates in metal-matrix composite materials; 8. The imaging of electronic device elements for future computing schemes, such as quantum computing. This application is particularly promising, since the ability to prepare the elements lithographically provides the a-priori information needed to solve the phase problem. In life sciences such a technique is needed to determine the bulk (internal) structure of assemblies of macromolecules (molecular machines), protein complexes, and virus particles at a resolution sufficient to recognize known proteins and determine their relationships to each other. The apparatus we have developed for CXDI is reviewed in section 3 and reconstruction methods are described in section 4. Section 5 summarizes data on the radiation damage limitations to resolution, and the theoretical basis for these limits, while Section 6 describes a method for overcoming those limits by using ultrafast x-ray pulses. We begin by describing a new algorithm for CXDI which avoids the need for a priori knowledge of the object support.

## II The CXDI technique

A CXDI experiment consists of three steps: (a) the sample is illuminated by monochromatic coherent x-rays and a recording is made of a single diffraction pattern (for 2D) or a tilt series (for 3D); (b) the phases of the pattern are recovered from the measured intensities using established phase-retrieval algorithms; (c) the unknown object is then recovered by Fourier inversion. In the Gerchberg-Saxton-Fienup scheme one starts with random phases, which lead to noise when transformed from reciprocal to real space. One then imposes the “finite support” constraint (namely that there must be a blank frame around the specimen), before transforming back to reciprocal space. In reciprocal space the phases so generated are combined with the measured diffraction magnitudes to start the next iteration. After a large number of iterations, in most cases the object emerges from the noise. Our novel improvement is that the estimate for the object support is continually updated by thresholding the intensity of the current object reconstruction marchesini:2003 . We start from a threshold of the transform of the diffraction pattern and as the iterations progress the support converges to a tight boundary around the object. This, in turn, improves the image reconstruction, which gives a better estimate of the support. An example of the reconstruction of a simulated diffraction pattern produced by a cluster of gold balls is shown in the movie (Fig. 3) together with the support.

The algorithm does not require any “atomicity” constraint provided by the gold balls as demonstrated by the reconstruction of a greyscale image. The algorithm also successfully reconstructs complex objects (those that cause large variations in the phase of the exit wavefield in two dimensions), which hitherto have been experimentally difficult to reconstruct. This opens up the possibility of image reconstruction from microdiffraction patterns, where the illumination is tightly focused on the object.

## III CXDI experiments at the ALS

Our experiments in coherent diffraction began with trials using electron weiestrall:2002  and visible-light spence:2002  optics and continued with CXDI at the Advanced Light Source (ALS). The experiments used the “pink” beam at beam-line 9.0.1 howells:2002  which is fed by a 10-cm-period undulator operating in third harmonic with deflection parameter ($`K`$) equal to 1.2 and delivering 588 eV (2.11 nm) photons. Features of the beam line include a 0.5 $`\mu `$m-thick, 750 $`\mu `$m-square Be window to separate the UHV beam line from the low-vacuum sample environment, a monochromator consisting of an off-axis segment of a zone plate and the diffraction experiment itself (Fig. 1).

The x-ray coherence length $`l_c`$ must be greater than the maximum path difference between any pair of interfering rays, i. e. $`l_c>w\theta _{\mathrm{max}}`$, where $`w`$ is the width of the sample and $`\theta _{\mathrm{max}}`$ is the maximum diffraction angle. For our geometry and wavelength, $`\theta _{\mathrm{max}}`$= 0.12 radian and the resolution limit is 8.4 nm. For the 5 $`\mu `$m aperture (effectively the monochromator exit slit) shown in Fig. 1, the resolving power is about 500, the coherence length is then 1$`\mu `$m and the maximum sample illumination area 8$`\times `$8 $`\mu `$m2. Similarly the (spatial) coherence patch provided by the 5 $`\mu `$m aperture is 10$`\times `$10 $`\mu `$m2. Allowing for an empty (but still coherently illuminated) band around the sample, its allowed size is thus $`<`$ 4$`\times `$4 $`\mu `$m2.

We consider now the sampling of the diffraction pattern. The Shannon interval for frequency-space sampling of the intensity is $`1\mathrm{/}\left(2w\right)=\mathrm{\Delta }\mathrm{/}\lambda z`$ where $`z`$ is the sample-to-detector distance and $`\mathrm{\Delta }`$ is the detector-plane increment (a 25 $`\mu `$m CCD pixel in our case). For our $`\lambda `$ and $`z`$ values this also leads to a maximum sample width of 4 $`\mu `$m. This is correct (Shannon) sampling of the diffraction-plane intensity and twofold oversampling in each direction of the diffraction-plane wave amplitude. Note that these limits on the sample size arising from coherence and sampling considerations are not the only ones in effect.

We have carried out three series of experiments, all using test samples made from 50 nm gold balls. The first he:acta:2003  demonstrated the basic 2D technique with image reconstruction using a support function determined by scanning electron microscopy. The second he:acta:2003  used a sample intentionally prepared in two separated parts, and reconstruction was achieved using information from the 2D diffraction pattern alone. The third series used a miniature sample-rotation device to collect several tomographic data sets. The picture of the device shown in Fig. 1 (picture width = 7 cm) shows the sample and its rotation spindle and driver. The small black square is the frame of a Si3N4 window on which the 3D sample is deposited in a 2.5-$`\mu `$m-wide microfabricated tetrahedron (Fig. 2). Not shown is an angular Vernier scale that was used to measure the rotation angle. Using this apparatus, a set of 150 views with at least a 100 second exposure time per view required about 10 hours. The 3D data generated by the object shown in Fig. 2 are still being analysed.

## IV Reconstructions

Samples were made by placing a droplet of solution containing ‘gold conjugate’ colloidal gold balls on a silicon nitride window (thickness 100 nm) and allowing it to dry. The gold balls formed several single layered (2D) clusters on the SiN membranes, as imaged by a field-emission scanning electron microscope (SEM).

In the first experiment the use of a Si3N4 window of the order of 5 microns width ensured that the sample was isolated and of the required size he:acta:2003 . The structure contained at least one isolated ball generating a ‘reference wave’ which interfered with the signals from other clusters to form a hologram.

The autocorrelation function obtained by Fourier transforming the intensity of the diffraction pattern included an image of every cluster convolved with the single ball, and these images formed a faithful representation of the structure in real space, without iterative processing. Not all the clusters could be imaged this way, since some of the intra-cluster distances were overlapping. The complex transmission function of the Si wedge at the corner of the window generated a complex object that was difficult to reconstruct by phase retrieval. For real objects, we have seen in simulations that using a support of triangular shape was sufficient to obtain the image. In the experimental case, we had to use a support obtained from a low resolution image of the object obtained by scanning electron microscopy (SEM). In the next experiment the object was isolated by means of sweeping the sample particles with an Atomic Force Microscope (AFM). The reconstruction has been obtained using the standard HIO algorithm fienup:1982  with missing data due to the beam stop and its supports left unconstrained during the iterations. Rather than rely on a low-resolution secondary imaging method to obtain an estimate of the support, we obtained an estimate of the object support from the support of its autocorrelation. The object support was then adjusted manually as the reconstructed image became clearer he:prb:2003 . A recently developed algorithm allows both initial support selection and adjustment to be performed automatically marchesini:2003 , extending the technique to objects other than a few separated clusters. This algorithm has been successfully applied on simulated 3D diffraction patterns, and it is being tested on the experimental 3D data recently recorded. We are also investigating the use of the real space algorithm hauriege:2003  to reconstruct the clusters in Fig. 3 and the 3D data. While this algorithm does require prior information in the form of a low-resolution image, it may be advantageous for the reconstruction from sparse, irregular, incomplete and noisy data.

## V Dose and Flux limitations to performance

In order to assess the promise of CXDI we have carried out calculations intended to determine the fluence (total photons per unit area) and dose (absorbed energy per unit mass) required to make a 3D CXDI image at given resolution. The basis of the calculation was the so-called “dose fractionation” theorem hegerl:1976 ; mcewen:1995 . The theorem states that the radiation dose required to measure the strength of a voxel is the same for a multi-view 3-D experiment (with appropriate reconstruction algorithms) as it is for measurement of the voxel alone, provided that the voxel size and the statistical accuracy are the same in both cases. The results are therefore based on a calculation of diffraction by a single voxel. The conclusions were as follows:

1. The dose ($`D`$) and fluence ($`N`$) needed to produce $`P`$ scattered x-rays per voxel into a detector collecting the angle required for resolution $`d`$, are given by $`N=P/\sigma _s`$ and $`D=\mu h\nu P/\epsilon \sigma _s`$ where $`\mu `$ is the absorption coefficient, $`h\nu `$ the photon energy and $`\epsilon `$ the density.
2. The scattering cross section $`\sigma _s`$ of the voxel of size $`d\times d\times d`$ is given by $`\sigma _s=r_e^2\lambda ^2\left|\rho \right|^2d^4`$ (see also henke:1955 ) where $`r_e`$ is the classical electron radius, $`\lambda `$ the wavelength and $`\rho `$ the complex electron density.
3. For resolution $`d`$, the dose and fluence scale as $`d^4`$
4. For light elements $`\rho `$ is fairly constant above about 1 keV so the energy dependence of the dose is expected to be quite flat in the range 1- 10 keV.
5. $`N`$ is dominated by the cross section and scales with the square of photon energy. Moreover, the coherent power of a source of brightness $`B`$ is $`B\left(\lambda /2\right)^2`$. This implies a fourth-power fluence penalty for increasing the x-ray energy.
6. Therefore one should use the lowest possible x-ray energy consistent with (roughly) $`\lambda <d/2`$
7. The dose for detecting a 10 nm protein feature against a background of water according to the Rose criterion rose:1948  is shown in Fig. 4. The required imaging dose in the energy range 1-10 keV is roughly $`10^9`$ Gy

Quantitative dose limit: to obtain an estimate of the resolution limit we have plotted a variety of literature values of the dose needed to destroy biological features as a function of feature size (Fig. 1). Also plotted is the required imaging dose. The tentative conclusion from the graph is that the resolution of the crossover, i.e. about 10 nm, should be possible for unlabelled life-science samples although for material-sciences samples the radiation tolerance (and thus the resolution) can be much higher.

Quantitative fluence limit: in our latest experiment at ALS, we collected a full 3D data set at a resolution that we believe to be around 10 nm (although this is not yet supported any reconstructions) in about 10 hours. We project that a beam line optimized for this experiment operating on the ALS after its planned performance upgrade would collect diffraction data about $`10^4`$ times faster than now. From 2, above, this should allow the step from 10 nm to around 1 nm resolution for sufficiently radiation-resistant samples.

## VI Femtosecond CXDI, beyond the radiation-damage limit

A way of overcoming the radiation damage limit in x-ray imaging is to use pulses of x-rays that are shorter in duration than the damage process itself. This idea of flash imaging, first suggested by Solem solem:1982 ; solem:1984 , has been proposed to be extended all the way to atomic resolution using femtosecond pulses from an x-ray free-electron laser (XFEL) neutze:2000 . The methodology of CXDI could be used in this case to image single molecules neutze:2000 ; szoke:1999 ; miao:2001 . The general concept for imaging non-periodic samples is to inject reproducible samples (macromolecules, complexes, or virus particles) into the XFEL beam, to intersect one particle (and record one diffraction pattern) per pulse.

In the general XFEL experiment the particle orientation in three dimensions will be random and unknown, and the individual diffraction patterns will be noisy (especially at the highest diffracted angles where the highest-resolution orientation information resides). General ideas and methods, developed for “single-particle” reconstructions in cryo-electron microscopy franck:1996 ; vanheel:2000 , for extracting tomographic information from a huge ensemble of randomly oriented noisy images can be applied here. Just as in single-particle EM, the limitations to image resolution are the ability to sort patterns into classes of like-oriented particles, so they can be averaged to improve signal to noise, and the reproducibility of particles. A statistical analysis huldt:2003  has shown that signal levels much less than a single count per pixel are adequate to be able to classify diffraction patterns. Even so, due to the exceedingly small scattering cross section of a single macromolecule, pulse fluences of $`>`$1012 photons per (100 nm)2 are required to achieve the required signal levels. This corresponds to a dose per particle of $`>10^{14}`$ Gray. The dose could be reduced substantially if methods could be employed to orient the particle in free space (even if only along a single axis) or if symmetric nanoclusters of particles could be formed.

We have performed calculations to determine, given the required fluence, the longest pulse duration possible to acquire the diffraction pattern before disruption of the electron density of the particle at the atomic scale. The calculation uses a hydrodynamic method treating the electrons and the ions as two separate fluids that interact through the Coulomb force and ionization processes, the assumption being that a large enough macromolecule can be treated as a continuum. Although it does not treat the atomic motions as accurately as the molecular dynamics model neutze:2000 , the hydrodynamic model is computationally faster and it also enables the inclusion of additional physics effects, in particular the inclusion of the redistribution of free electrons throughout the molecule.

The dominant interactions between the matter and x-ray pulse are $`K`$-shell photoionization, (producing damage) and elastic scattering producing the diffraction pattern. Following x-ray absorption in about 2 to 10 fs mcguire:1969 , the $`K`$-shell holes decay by emitting Auger electrons with energies in the range of 250 to 500 eV furguson:1989 . As the photo and Auger electrons escape they leave behind a molecule with increasing charge and the growing population of low energy secondary electrons becomes trapped and grows still further by collisional ionization (causing the bulk of the damage).

The trapped electrons quickly relax in energy and position to form a neutralizing cloud around the positively charged ions. The particle assumes roughly a two-zone structure, consisting of a neutral core and a positively charged outer shell similar to Debye shielding. On a longer timescale of several 10’s fs, a macroscopic motion of the whole molecule, called a Coulomb explosion, takes place, leading to an outward radial motion of the ions. The ion motion is greatest within the charged outer layer of the molecule. In the core of the molecule the electron temperature is highest, with the greatest amount of ionization and blurring of the electron density. Preliminary results from the hydrodynamic model indicate that collisional ionization is very rapid and may limit the maximum pulse length to a value smaller than the 10-20 fs inferred by considering only photoionization, Auger ionization and atomic motion. The limits imposed by collisional ionization might be overcome by developing a method to reconstruct atomic positions from partially ionized atoms.

To date it has not been possible to acquire experimental data on the rate of the Coulomb explosion and the effect on electron density of ultrafast high-intensity x-ray pulses. We plan to investigate this at the Tesla Test Facility Phase 2 tesla:2003 , a soft-x-ray FEL that will be operational in 2004 initially at a wavelength of 20 nm. The first experiments that are planned are to determine the onset of the Coulomb explosion of particles in an aerosol, as a function of pulse fluence, as measured by diffraction from the particles by the x-ray pulses themselves.

## VII Conclusions

We have discussed advantages and limitations of the CXDI technique, and our progress in implementing the experiments and reconstruction techniques. Preparation of special objects containing an isolated gold ball near an unknown compact object are shown experimentally to allow simple (although low resolution) image reconstruction without iteration based on the autocorrelation function (Fourier transform holography). Such an image may also provide a support function for the retrieval of higher-resolution images using the HIO algorithm. We have developed an algorithm that obtains the support automatically from the autocorrelation function by directly Fourier transforming the recorded diffraction intensity. This eliminates the need for an auxiliary imaging experiment to obtain the support function.

Our last experiment shows that the exposure times for 10 nm resolution are reasonable (about 10 hours) and this resolution is near the dose limit for life-science experiments. Thus we conclude that 3D 10-nm resolution life science experiments on frozen hydrated samples will be possible, and even better resolutions might be achievable with labeling techniques. Calculations indicate that a factor of about $`10^4`$ improvement in imaging speed can be expected using a dedicated CXDI beam line at the upgraded ALS. This is important for applications in material science where the dose limitation is much less severe. Based on the inverse fourth power scaling of the required flux with resolution, this should allow a factor ten improvement in resolution to about 1 nm.

These results have prompted us to develop two parallel projects. One is the development of a program for cryo-CXDI at ALS, and the other aims at femtosecond imaging of single molecules at the LCLS and other $`4^{th}`$ generation sources.

The work in CXDI at ALS will be performed in collaboration with groups from NSLS and Stony Brook, who have constructed an instrument designed for collecting 3D data sets beetz:2003 , and which will be installed at ALS. The instrument incorporates a set of in-vacuum stepper motors for precision positioning the required apertures, beam stop, and diagnostics. It incorporates an air-lock and precision goniometer (developed for electron tomography) for easy introduction and manipulation of the specimen. The specimen holder is designed so as to be able to rotate the specimen over the angular range from -80 to +80 degrees before the specimen supporting grid obscures the beam. A unique feature of the instrument is that a zone plate can be positioned between the specimen and the CCD, thereby a low resolution image of the specimen can be recorded. This feature helps align the specimen, and provides information on its “support”.

The key idea to achieving atomic resolution imaging of radiation-sensitive single molecules and particles is that the damage limit may be overcome with the use of very short x-ray pulses to capture the data before damage occurs. Three problems need to be assessed to be able to perform single molecule imaging experiments: femtosecond damage and experiment modeling; image orientation, characterization and reconstruction; and single molecule sample handling.

Fast and accurate hydrodynamic models to describe the interaction of a molecule with femtosecond x-ray pulses needs confirmation from experimental results of short pulse, high-field x-ray-matter interaction experiments planned in the near future. The recorded diffraction patterns must be classified, averaged to increase the signal-to-noise ratio, and oriented for the final 3D reconstruction. Methods developed for single-particle electron microscopy can be applied here. We must also develop methods to inject samples into the beam and, if possible, orient the particles (or at least influence their orientation). Preliminary lower resolution experiments on single particles will provide the answers to the experimental needs to build an experiment at LCLS to perform single-molecule imaging.

We have established the fundamental boundaries of applicability of CXDI to problems in materials and biological sciences. With the rapidly growing importance of nanostructures, and the potential that future developments in this area have for major breakthroughs in the fundamental and applied sciences and in technology, the addition of a new probe that can image thick objects at very high 3d spatial resolution will have a decisive impact on nanoscience and technology. Even more ambitiously, atomic-resolution CXDI at an XFEL, has the potential to provide atomic structure determination of biological systems.

###### Acknowledgements.

We would like to thank Rick Levesque of LLNL for designing and building the rotation stage, and Cindy Larson and Sherry Baker (both LLNL) for SEM imaging. This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under Contract No. W-7405-ENG-48 and the Director, Office of Energy Research, Office of Basics Energy Sciences, Materials Sciences Division of the U. S. Department of Energy, under Contract No. DE-AC03-76SF00098. SM acknowledges funding from the National Science Foundation. The Center for Biophotonics, an NSF Science and Technology Center, is managed by the University of California, Davis, under Cooperative Agreement No. PHY0120999.


text:
# THE FORMATION AND EVOLUTION OF MASSIVE STAR CLUSTERS: HISTORICAL OVERVIEW

## 1 INTRODUCTION

### 1.1 Why there are star clusters and associations?

The existence of multiple stars and star clusters is, as Larson (2003) has recently reminded us, due to the fact that the angular momentum in typical dense cores of giant interstellar molecular clouds (GMCs) is three orders of magnitude greater than can be contained within a single star. The magnetic fields in GMCs cannot dispose of all of this angular momentum. This is so because it has become both dynamically unimportant, and decoupled from the gas, during the later stages of the collapse of these GMC cores. The dense cores of such molecular clouds are therefore expected to fragment into clumps of massive stars. Such concentrations of young stars were named “associations” by Ambartsumian (1954). Associations are now known to come in three flavors: (1) OB-associations, which are comprised of massive young main sequence stars, (2) T-associations, which consist of young intermediate-mass stars that are still contracting to the main sequence, and (3) R-associations (van den Bergh 1966), which are made up of stars that illuminate the dusty interstellar clouds in which they are still embedded. Since the space density of R-associations is much greater than that of OB-associations such R-associatiations are particularly suitable for the mapping of nearby Galactic spiral structure (van den Bergh 1968).

### 1.2 Structure of associations.

Within low-density expanding (positive energy) associations one often also also finds denser gravitationally stable (negative energy) star clusters (e.g. Blaauw 1964). Some of these clusters, like the small (R = 0.3 pc) clustering near o Persei (that is located in the emission nebula IC 348) are so small that they clearly have not yet expanded out to their eventual equilibrium radii. Not all associations have the same structure. The 30 Doradus association (see Fig. 1 of Walborn, Maíz-Apellániz & Barbá 2002) in the Large Magellanic Cloud is centered on the single massive compact young OB cluster R136 (which forms the core of the 30 Doradus association), whereas the comparably rich association NGC 604 in M33 (see Fig. 2 of Miskey & Bruhweiler 2003) contains many smaller subclusterings, rather than a single massive core. NGC 206, which is the most luminous association in the Andromeda galaxy (M31), is a slightly older example of a large association that does not have a single massive core. It is presently not known why some associations exhibit such subclustering whereas others do not.

## 2 STAR CLUSTERS AND GALACTIC EVOLUTION

### 2.1 Ages of halo and disk clusters.

Hansen et al. (2002) have found that the faintest white dwarfs in the globular cluster M 4 are 2.5 mag fainter than are the faintest field white dwarfs in the Solar neighborhood. The resulting cooling ages are 12.7 $`\pm `$ 0.7 Gyr for M 4 and 7.3 $`\pm `$ 1.5 Gyr for the Galactic disk. The relatively young age of the nearby disk has recently been confirmed by Sandage et al. (2003) who determined the age of old field subgiants from their Hipparcos parallaxes. These results show, beyond reasonable doubt, that the Galactic disk is $``$5 Gyr younger than the Galactic halo. As was first emphasized by Berman & Suchkov (1991) the only reasonable way to obtain such a large age difference is by stopping star formation via heating of the Galactic gas to temperatures in excess of that defining virial equilibrium. Such heating would likely have resulted from the strong stellar wind generated by the burst of star formation that accompanied the formation of the Galactic bulge. The fact that the oldest open clusters, such as NGC 6791 and NGC 188, have metallicities that are 1 - 2 times Solar shows that these objects must have been formed from gas that had already been greatly enriched by star and supernova formation.

### 2.2 Cluster formation and galactic star formation rates.

Larsen & Richtler (1999, 2000) have recently shown that the fraction of the total light of galaxies that is locked up in star clusters differs widely from galaxy to galaxy. Star clusters are, for example, almost absent from the very inactive Local Group galaxy IC 1613 (Baade 1963). On the other hand 5.6% of the total V-band light of the nearby starburst galaxy NGC 1569 is emitted by stars in clusters. If massive clusters are mainly formed during intense star bursts then : (1) The history of cluster formation in a galaxy cannot be used as a proxy for the overall history of star formation in that galaxy, and (2) galaxies that exhibit an above average specific cluster frequency (Harris & van den Bergh 1981) probably had at least one major starburst in their history. In view of the apparent correlation between starbursts and the formation of massive clusters it would be particularly interesting to see if the specific cluster frequency is enhanced in IC 10, which is the nearest starburst galaxy. It seems likely (Elmegreen & Efremov 1997, Tan & McKee 2002) that massive clusters of all ages preferentially formed (and survived) in high pressure environments. This is so because virialized clouds are more tightly bound at high pressures than they are in low pressure environments. Such high pressure might be due to either a high background virial density or to large-scale shocks that are generated during galaxy collisions. It would be interesting to know if the great burst of globular cluster formation, that took place $``$12 Gyr ago, was also triggered by collisions, or if shocks associated with the reionization of the Universe (See Miralda-Escudé 2003 for a review) might have triggered the great early burst of globular cluster formation.

### 2.3 Cluster formation and galaxy interactions.

Van den Bergh (1960) found that the spiral arms of some late-type interacting galaxies were characterized by particularly patchy (DDO type Sc\*) spiral arms. Since such patches are due to the presence of numerous star clusters and associations, this indicated that strong tidal interactions and mergers can trigger bursts of cluster formation. The notion that interactions between galaxies might enhance the formation rate of massive star clusters (e.g. Schweitzer 1987) has lately received strong support from observations with the Hubble Space Telescope \[see Whitmore (2003) for a recent review\]. The fact that the fraction of all star formation that takes place in massive clusters is apparently enhanced during bursts of star formation indicates that the history of cluster formation in a particular galaxy might be quite different from the overall history of star formation in that same galaxy. In particular it appears likely that the specific cluster formation frequency was enhanced during the early phase of galaxy formation that was dominated by mergers. Proto-planetary disks surrounding young stars might have a high probability of being disrupted in cluster stars than in field stars. It is therefore quite possible that the fraction of stars that form planets may be reduced in early generations of star formation. However, it will be difficult to test this prediction because the low metallicity of the first generation of stars might also inhibit planet formation. A possible counterexample is provided by the old planet that Sigurdsson et al. (2003) have recently found in the globular cluster M 4. Nevertheless the fact that the disk fraction in very young clusters decreases with age (Lada & Lada 2003) does show that the cluster environment can destroy protoplanetary disks on a time scale of a few million years.

### 2.4 Formation of star clusters and cluster location.

Morgan (1958, 1959) devised a new one-dimensional system of galaxy classification, in which the principal classification parameter was central concentration of light in the galaxy image. In the footnotes to his main tables Morgan drew attention to the fact that many spiral galaxies have knots of very active star formation just outside their nuclei. The Arches (age 2 Myr) and Quintuplet (age 4 Myr) are examples of such recently formed super star clusters in own Milky Way system. It is presently not clear why the formation of massive star clusters is favored in a region where differential rotation and tidal stresses would appear to make the formation (and subsequent gravitational collapse) of the massive cores giant molecular clouds particularly challenging. Barth et al. (1995) find that very luminous star clusters may also occur at somewhat greater nuclear distances in the star forming rings that surround some active active galactic nuclei.

### 2.5 Luminosity functions of stars in young star clusters.

Larson (1999) has recently reviewed the observational evidence on the luminosity function of individual stars in young star clusters. He concludes that “No clear evidence has been found for any systematic dependence of the IMF \[initial mass function\] on any property of the system studied, and this has lead to the current widely held view that the IMF is universal, at least in the local universe.” However, the Arches cluster near the Galactic center (Figer et al. 1999) may have a flatter mass spectrum than typical young clusters in the Solar neighborhood. So, perhaps, the mass spectrum of star formation in clusters is not entirely independent of environment after all.

### 2.6 The mass spectrum of cluster formation.

In his book “Sternhaufen”, which is the first monograph entirely devoted to star clusters, ten Bruggencate (1927) divides clusters into three classes: (1) globular clusters, (2) open clusters, and (3) star swarms, a class that would, in modern terminology, consist of both moving groups and stellar associations. Ten Bruggencate emphasized the fact that globular clusters have a spheroidal distribution that is centered on the Galactic bulge, whereas open clusters are mainly located in the flattened disk of the Milky Way. Furthermore he pointed to the differences between the color- magnitude diagrams of these two types of clusters, which we now know to be due to systematic differences in age and chemical composition. Perhaps surprisingly, ten Bruggencate did not discuss the systematic difference between the integrated luminosities of typical globular and galactic clusters. Modern work broadly supports the division of Galactic star clusters into (1) old, luminous, metal-poor globular clusters and (2) young, less luminous, and metal-rich open clusters. However, the more extensive data base that is presently available does show intermediate objects such as (a) metal-rich globulars in the Galactic bulge, (b) a few open clusters such as NGC 6791 and NGC 188 with ages of 5-10 Gyr, and (c) some outer halo globular clusters with an ages on only $``$10 Gyr. Surprisingly the age dichotomy between globulars and open clusters is even more pronounced in the Large Magellanic Cloud (Da Costa 2002) than it is in the Galaxy. On the other hand the systematic difference between the masses (luminosities) of open and globular clusters is less pronounced in the Large Cloud than it is in the Milky Way, i. e. the LMC open cluster population contains a larger fraction of massive clusters than does that of Galactic open clusters.

## 3 THE EVOLUTION OF CLUSTER SYSTEMS

### 3.1 Unimodal and bimodal metallicity cluster systems.

It has been known for may years (van den Bergh 1975) that the mean metallicities of globular cluster systems grows with increasing parent galaxy luminosity. This shows that the nature of globular cluster systems is inextricably linked to the characteristics of their host galaxies. This leads one to suspect that the properties of globular cluster systems might also correlate with the Hubble types of their parent galaxies. Unfortunately the number of globular cluster systems for which detailed information on the metallicities of individual clusters has become available is still small. The data are particularly incomplete for spiral galaxies. However, Eerik & Tenjes (2003) have recently published quite extensive information on the metallicity distributions of globular cluster systems surrounding E7-E7 and S0 galaxies. Their data appear to show (97.5% confidence) that elliptical galaxies are more likely to have bimodal cluster metallicity distributions than is the case for S0 galaxies. However, this effect appears to be due to the fact that that S0 galaxies are systematically fainter than ellipticals. On average ellipticals therefore have richer globular cluster systems than do lenticulars. The observed effect may be due to the fact that it is difficult to establish bimodality of the cluster metallicity distribution in a poor cluster systems. It would be very interesting to extend the study of the relationship between modality of cluster systems to spiral galaxies. A practical problem is, however, that spirals (1) have a lower specific globular clusters frequency than do ellipticals, and (2) that it is often difficult to distinguish unambiguously between stars and globular clusters in distant spirals.

### 3.2 The flattening of globular clusters.

The majority of both open and globular clusters in the Galaxy appear as little elongated almost spherical objects. A strikingly different situation is observed to hold in the Large Magellanic Cloud (Frenk & Fall 1982, van den Bergh & Morbey 1984) in which both globular clusters and open clusters are found to be systematically more flattened than are those in the Milky Way System. Such highly flattened objects as NGC 121 in the SMC and NGC 1978 in the LMC are very rare or absent in the Galaxy. The reason for this systematic difference is presently not known. In this connection it is also of interest to note (van den Bergh 1996) that the most luminous globular in a cluster system (e. g. Omega Centauri in the Galaxy, Mayall II in M31, NGC 1835 in the LMC and NGC 121 in the SMC) is often also the most highly flattened.

## 4 DESIDERATA FOR FUTURE WORK

The topics covered in this review suggest that it might be rewarding to study the following problems in more detail:

* Why is the age distribution of LMC globular clusters so different from that of Galactic globulars? This difference is particularly puzzling because the dichotomy between the masses of open and globular star clusters appears to be much less pronounced in the Large Cloud than it is in the Galaxy.
* We need to understand why the evolutionary histories of different globular cluster systems differ. Why do some have bimodal metallicity distributions, whereas others do not.
* Why are there systematic galaxy-to-galaxy differences between the average flattenings of individual star clusters? Why do these mean flattenings depend on cluster luminosity, but not on cluster age?
* Some evidence suggests that the mass spectrum with which stars form in clusters is universal, but other data indicate that some clusters may have been formed with “top heavy” mass spectra. Resolving this problem would help us understand both the dynamical and chemical evolution of star clusters.

I also belatedly wish to thank Prof. ten Bruggencate for giving me his personal copy of his book “Sternhaufen” just before his death. I am also indebted to Doug Johnstone for helpful discussions.