text:
# Evidence for parametric memory and correlation in semiconductor microcavities

## Abstract

We measured the dynamics of polariton parametric stimulated scattering in semiconductor microcavities, by time-resolving the amplified signal with 250 fs resolution. Our experiments demonstrate that the stimulation process is considerably delayed with respect to the arrival time of both probe and pump pulses. This effect is clearly observable in our high quality sample due to the long lifetime of the microcavity polaritons (about 10 ps), and it is in excellent agreement with the model describing the coherent conversion of two pump polaritons into a signal-idler pair. We demonstrate that the non-instantaneous character of the polariton stimulation occurs because the polariton-wave amplification requires the build up of correlated signal-idler polariton fields.

Matter-wave amplification is a fascinating subject, which is attracting great interest in the physics of atomic Bose-Einstein condensates DengFWM ; InouyeFWM ; Vogels , but also in the field of semiconductor quantum optics Yamamotoboser ; Lesidang ; Rubo ; Messin . In semiconductor microcavities, the particles involved in the coherent collisions are exotic bosons, the so-called microcavity polaritons, resulting from the strong exciton-photon coupling Weisbuch . Recently, spectacular amplification phenomena due to stimulated polariton scattering have been experimentally demonstrated SavvidisPPA ; Sabanature ; Stevenson ; DasbachGain . Polariton-wave collisions can also be manipulated by judiciously phase-locked excitation pulses, which control the emission over the whole momentum space KundermannCC .

Particularly interesting is the parametric nature of the scattering. The basic underlying process is the coherent conversion of two pump polaritons into a signal-idler polariton pair conserving total energy and in-plane momentum thanks to the very peculiar shape of the dispersion relation of polaritons CiutiOPA . Following this description, the signal idler polariton pairs should be entangled Schwendimannstatistics ; SavastaPRL . Squeezing of the polaritons has been recently been observed in a geometry where signal and idler are degenerate which confirms the possibility to achieve quantum correlation by this non-linear effect Baas . The process is initially started by providing a small incident signal occupation with a weak resonant probe laser pulse. The dynamical aspects of the parametric scattering process are therefore particularly interesting, and very seldomly studied Erland2 . In contrast to the stimulation processes which occur for example in a laser, where the scattering rate into the one-particle final state is directly proportional to its occupation, the situation in the polariton system is more complicated. The stimulation of the parametric scattering involves the entire final states of the conversion process, i. e. the signal-idler pairs Schwendimannstatistics . The behavior should differ largely from the dynamics observed in the case of a laser (stimulation from an incoherent reservoir) and be similar to what occurs in a conventional optical parametric oscillator (coherent wave mixing). However, in our microcavity the parametric process involves real mixed states, whose excitonic part also influences the dynamics.

The investigation of the polariton stimulated dynamics is a challenging task, because the polariton signal temporal dynamics needs to be measured with a time-resolution much shorter than the relaxation times of the polariton quasi-particles. In order to achieve this goal, we have implemented an up-conversion set-up with sub-picosecond resolution ($`0.25`$ ps) to directly measure the temporal dynamics of polariton signal amplification. The experiments have been performed on a very high quality microcavity sample with very slow decoherence, whose polariton radiative lifetime exceeds $`10`$ ps. To our knowledge, our experiments demonstrate very directly for the first time the non-instantaneous nature of the polariton parametric stimulated scattering. Our claim is further susbstantiated by an excellent agreement with the predictions given by the polariton parametric amplifier model CiutiOPA  and correlation effects between signal and idler are evidenced.

The sample is a single quantum well $`\lambda `$-cavity held in a helium-bath cryostat at a temperature of 2 K. The quality of the sample can be assessed by either the linewidth of the normal modes, or the decay of the lower polariton polarization after pulsed excitation. The linewith is of the order of $`100\mu eV`$ and the photon lifetime in excess of 10 ps. The set-up is an angle-resolved pump and probe apparatus combined with an up-conversion detection (see Fig.1). The probe beam impinges onto the sample at normal incidence, whereas the pump arrives at an off-normal angle of about 10 deg., which is the so-called magic angle for producing parametric scattering for this microcavity. Both pump and probe are pulsed with 80 MHz repetition rate; the pump is spectrally filtered to selectively excite the lower polariton branch ($``$1 meV spectral width and $``$1 ps temporal duration), while the weak probe keeps a broader spectrum ($``$ 15 meV) and shorter duration ($``$150 fs). For the upconversion detection, the transmitted probe signal is focused onto a BBO crystal together with a 150 fs non collinear gating pulse. The cross-correlation beam generated in the crystal is then selected by a diaphragm and detected with a bi-alkali photomultiplier, so that the shape of the probe transmission in real time is reconstructed with a resolution around 250 fs by scanning the delay of the gate pulse. All experiments are carried out far from the regime of saturation, and clearly in the strong coupling regime.

The upper panel of Fig. 2 shows the measured temporal evolution of the stimulated parametric scattering for different pump densities. The origin of the time scale corresponds to the arrival time of the pump lasernoise . The probe hits the sample 1.5 ps before the pump (similar results are obtained for other delays between the pump and the probe). Curve (a) shows the temporal evolution of the transmitted probe signal in absence of the pump. Due to the simultaneous excitation of the upper and lower polariton branch the temporal evolution exhibits oscillations due to quantum beats between lower and upper polariton. This signal decays with the expected cavity lifetime. When the pump density is adjusted above threshold (trace 2 (b)), the probe signal first starts decaying before a rise due to the stimulated scattering amplification takes place. This rise is significantly delayed with respect to the arrival of the pump pulse. When the pump intensity is further risen, the scattering even saturates (2c, 2d). As previously noticed the process is significantly delayed with respect to the arrival of the pump.

Curves (b)-(d) in Fig. 2 show the delayed build up of the scattered signal. To get a better insight into this behavior, the experiments were performed for different pump to probe delays as shown on the lower panel. The pump intensity is $`60W/cm^2`$, as for trace (c). For an easier reading of the data, the weak signal arising from the ”spontaneous” parametric luminescence of the pump measured with the pump alone (trace 2 (e)) has been subtracted in each case. In all cases, the stimulated emission is strongly delayed with respect to the incident probe pulse. When the probe hits the sample first ($`5.5ps`$ (j) and $`1.5ps`$ (i)) the transmitted signal first decays as expected and then gets amplified, but only with some delay with respect to the arrival of the pump. Curve (i) peaks at an earlier time than (j) because the polariton intensity at $`t=0`$ and $`k=0`$ is higher. Even when the probe hits the sample after the pump, the stimulated emission is clearly delayed (traces (f), (g), and (h)). These non trivial dynamics are somewhat astonishing as they cannot be explained using a simple exponential dependence of the stimulation. This demonstrates that the often assumed stimulated process following (N+1) where N is the sole population of the k=0 polariton state is far too simple. To understand in detail the dynamical behavior we now turn to some theoretical considerations.

The basic process for this kind of polariton amplification is considered to be the parametric conversion of two pump polaritons with wave-vector $`k_p`$ into a signal-idler pair, with conservation of both total energy and momentum. In its simplest version, the scattering in the system can be described by three modes only, at $`k=0`$, $`k=k_p`$ and $`k=2k_p`$ in a factorized mean-field approximation meanfield ; CiutiOPA . Since the parametric scattering involves states on the lower polariton dispersion, it is sufficient to look at the lower polariton branch only to understand the basic underlying dynamics. The use of proper polariton states for the lower branch allows then to write the dynamics of the system as a set of three coupled differential equations.

$`i\mathrm{}{\displaystyle \frac{P_0}{t}}=(\stackrel{~}{E}_0i\gamma _0)P_0+E_{int}P_{2k_p}^{}P_{k_p}^2+F_0(t)`$ (1)
$`i\mathrm{}{\displaystyle \frac{P_{k_p}}{t}}=(\stackrel{~}{E}_{k_p}i\gamma _{k_p})P_{k_p}+2E_{int}P_{k_p}^{}P_0P_{2k_p}+F_{k_p}(t)`$ (2)
$`i\mathrm{}{\displaystyle \frac{P_{2k_p}}{t}}=(\stackrel{~}{E}_{2k_p}i\gamma _{2k_p})P_{2k_p}+E_{int}P_0^{}P_{k_p}^2`$ (3)
where $`P_0`$, $`P_{k_p}`$, $`P_{2k_p}`$ are the mean lower polariton fields for signal, pump, and idler wave-vectors respectively and $`E_{int}`$ is the coupling energy due to polariton-polariton scattering potential. The quantities $`P_{0,k_p,2k_p}`$ are the rescaled polarizations $`P_k(t)=\frac{\lambda _X}{\sqrt{A}}p_k(t)`$ ($`p_k`$ are the polarizations). On the rhs of each of the three equations, the first term is the ’free term’, i.e. the evolution of the polarization of non-interacting polaritons having energy $`\stackrel{~}{E}(k)`$. The homogeneous k-dependent linewidth $`\gamma _k`$ determines an exponential decay rate of the signal ($`\gamma _k`$ is determined by the losses of polaritons through the cavity mirrors and also includes non-radiative losses). Note that the energy $`\stackrel{~}{E}(k)`$ includes a renormalization due to the polariton-polariton interaction and therefore is slightly higher than the energy $`E(k)`$ calculated in the absence of pump polaritons. The second term is the parametric scattering rate: for the signal (idler) it is proportional to the square of the pump polarization (i.e., proportional to the pump intensity) and to the conjugate of the idler (signal) polarization. The symmetry between signal and idler equations reflects the fact that for each polariton scattering from the pump down to the signal there is another one scattering up to the idler and vice versa. This is the origin of the correlation between signal and idler polaritons. The last term in the first two equations ($`F_{0,k_p}(t)`$) corresponds to the external driving electric field, i.e., the amplitude of pump and probe laser pulses.

A crucial physical consequence of parametric scattering is that the rate of scattering from the pump into the signal mode is proportional to the idler field. More directly, we can write

$$\left(\frac{P_0}{t}\right)_{par}=\frac{E_{int}}{i\mathrm{}}P_{2k_p}^{}(t)P_{k_p}^2(t),$$
(4)
The dynamical equation for the idler polariton field can be formally solved, giving the following analytical result :

$$P_{2k_p}(t)=\frac{E_{int}}{i\mathrm{}}_{\mathrm{}}^t𝑑t^{}P_0^{}(t)P_{k_p}^2(t^{})e^{\frac{i}{\mathrm{}}(E_{2k_p}^{LP}i\gamma _{2k_p})(tt^{})},$$
(5)
where $`E_{2k_p}^{LP}`$ is the idler polariton energy and $`\gamma _{2k_p}`$ the corresponding broadening. By inserting Eq. (3) into Eq. (4), we can get the expression for the signal parametric scattering rate in terms of the signal and pump polariton fields only, namely :

$$\left(\frac{P_0}{t}\right)_{par}=\left(\frac{E_{int}}{i\mathrm{}}\right)^2_{\mathrm{}}^t𝑑t^{}K_{par}(t,t^{})P_0(t^{}),$$
(6)
where the parametric memory kernel reads

$$K_{par}(t,t^{})=P_{k_p}^2(t^{})e^{+\frac{i}{\mathrm{}}(E_{2k_p}^{LP}+i\gamma _{2k_p})(tt^{})}P_{k_p}^2(t).$$
(7)
Eq. (6) shows that the rate of pump scattering into the signal mode is proportional to the signal field itself, i.e. the scattering is stimulated. However, the stimulation is not instantaneous, because the stimulated scattering rate depends on the value of the signal field at all times $`t^{}<t`$. This memory effect takes place, because the parametric stimulation requires the coherent and correlated build-up of the signal and idler fields. Instantaneous stimulation is recovered only when the parametric memory kernel kernelref  $`K_{par}(t,t^{})\delta (tt^{})`$. From Eq. (7), we understand that in order to experimentally detect the parametric memory effect, the time-resolution of the measurement must be much shorter than the decoherence times of the polariton fields.

In order to compare the measured dynamics, we have solved numerically the nonlinear equations for the polariton parametric amplifier. Since the short probe pulses also excite the upper polariton branch, we have extended equations 1-3. The calculations have been performed in the exciton-cavity photon basis, so that the calculation accounts for the lower and the upper polariton branch. As previously shown for the analytic derivation of the memory kernel (in the lower polariton basis), the same behavior can be found again in the numerical equations. First the external probe laser field creates a polariton occupation in the signal state. This polariton occupation together with the pump polaritons result in a build up of the idler polarization. Since the signal and idler equations are symmetric, the established idler polarization stimulates the pump polaritons to scatter to the signal state and the signal build up is delayed and non mono-exponential.

Fig. 3 shows the simulation of the dynamics using the model of the parametric amplifier as discussed before. The pump and probe densities have been adjusted according to the external pump and probe power by using a global conversion factor and the exciton and cavity linewidth have been inserted according to previous works on similar samples Stanleyquench . The pump to probe delays are the same as in the experiment (Fig. 2). The agreement with the experimental data is very good. All the features observed in the experiment (Fig. 2) are reproduced, even the saturation at the highest pump level. A very important observation is made looking at curves 2 (h) and 3 (f-h): the oscillations due to the LP-UP beating remain even after the pump has hit at t=0, which demonstrates that we are in the strong coupling regime in all cases. However the beat oscillations persist much longer in the simulation than in the experiment because of the absence of additional dephasing mechanisms for the upper polariton in the model dephasingmechanism . The saturation of the signal intensity observed for curves 2 (b, d) and 3 (b, d) is due to the depletion of the pump polariton reservoir because the scattering rate out of the pump reservoir increases with the square of the pump intensity (see equation 7) KundermannOECS03 .

The simulation also gives access to the pump and idler temporal evolution. To further inspect the dynamical interplay of signal, pump, and idler we look at the evolution of the three modes. In this simulation, which is depicted in Fig. 4, the probe hits the cavity again 1.5 ps before the pump.

The idler starts to build up when the pump hits and the signal follows the evolution of the idler after its initial decay. The simultaneous build up of signal and idler demonstrates their mutual correlation. This shows that even in the mean-field approximation a signal idler correlation appears in spite of the fact that all the field operators have been factorized. In Fig. 2 (e) we observe that the signal due to the ”spontaneous” parametric luminescence of the pump needs much more time to build up. The ”seed” for this signal is not the initial probe polariton occupation but the vacuum fluctuations of signal and idler modes (a model accounting for that can be found in Schwendimannstatistics ). Since the effect of the fluctuations starts to build up the signal much later (see curve 2 (e)), their effect can be neglected for the pump probe experiment, in which the correlations are built up very quickly due to the presence of the probe polaritons. The quantum beats between upper and lower polariton resonance appear very clearly in the signal, they also give rise to faint oscilltations in the pump and idler branches.

In conclusion, we have been able to observe the non-instantaneous nature of the polariton parametric scattering in semiconductor microcavities by time resolved pump-probe experiments. The excellent quality of the sample together with the high time resolution of the experimental setup allow to report this effect for the first time. The measurements are in excellent quantitative agreement with a mean field theory and it is shown that the stimulation takes place due to the mutual correlation of signal and idler. The demonstrated signal idler correlation makes the polaritons a promising system for emission of entangled photons and for quantum information processing.

We thank U. Oesterle, R. Stanley, R. Houdré, and M. Illegems who provided us with the excellent microcavity sample. Furthermore we would like to acknowledge fruitful discussions with A. Baas, R. Idrissi, A. Quattropani, V. Savona, and P. Schwendimann. The present work has been supported by the Swiss National Science Foundation.

Present address: Center for Ultracold Atoms, MIT, Cambridge, Ma, USA.


format: text

text:
# Qualms concerning Tsallis’ Use of the Maximum Entropy Formalism

## Abstract

Tsallis’ ‘statistical thermodynamic’ formulation of the nonadditive entropy of degree-$`\alpha `$ is neither correct nor self-consistent.

It is well known that the maximum entropy formalism , the minimum discrimination information , and Gauss’ principle  all lead to the same results when a certain condition on the prior probability distribution is imposed . All these methods lead to the same form of the posterior probability distribution; namely, the exponential family of distributions.

Tsallis and collaborators  have tried to adapt the maximum entropy formalism that uses the Shannon entropy to one that uses a nonadditive entropy of degree-$`\alpha `$. In order to come out with analytic expressions for the probabilities that maximize the nonadditive entropy they found it necessary to use ‘escort probabilities’ of the same power as the nonadditive entropy.

If the procedure they use is correct then it follows that Gauss’ principle should give the same optimum probabilities. Yet, we will find that the Tsallis result requires that the prior probability distribution be given by the same unphysical condition as the maximum entropy formalism and, what is worse, the potential of the error law be required to vanish. The potential of the error law is what information theory refers to as the error ; that is, the difference between the inaccuracy and the entropy. Unless the ‘true’ probability distribution, $`P=(p(x_1),p(x_2)\mathrm{},p(x_m))`$ coincides with the estimated probability distribution, $`Q=(q(x_1),q(x_2),\mathrm{}q(x_m))`$, the error does not vanish. Moreover, we shall show that two procedures of averaging, one using the escort probabilities explicitly, do not give the same result, and the relation between the potential of the error law and the nonadditive entropy requires the latter to vanish when the former vanishes.

Let $`X`$ be a random variable whose values $`x_1,x_2,\mathrm{},x_m`$ are obtained at $`m`$ independent trials. Prior to the observations the distribution is $`Q`$, and after the observations the unknown probability distribution is $`P`$. The observer has at his disposal the statistic

$$\widehat{a}=\frac{1}{m}\underset{i=1}{\overset{m}{}}x_i$$

to help him formulate a guess as to the form of $`Q`$. Gauss’ principle assumes that the probability distribution $`P`$ depends on a parameter $`a`$

$$a=E(X)=\underset{i=1}{\overset{m}{}}x_ip(x_i),$$
(1)
such that the arithmetic mean, $`\widehat{a}`$, is the maximum likelihood estimate of $`a`$. Furthermore, $`P`$ will depend upon the parameter $`a`$ in such a way that there is a value $`a^0`$ for which $`p(x_i;a^0)=q(x_i)`$, the prior distribution.

The maximum likelihood estimate,

$$\frac{}{a}\mathrm{log}(a)=0,$$

will lead to the exponential family of distributions when the log-likelihood function

$$\mathrm{log}(a)=\underset{i=1}{\overset{m}{}}\mathrm{log}p(x_i;a).$$

The likelihood equation

$$\frac{}{a}\psi (x_i;a)=0,$$

where $`\psi (x_i;a)=\mathrm{log}p(x_i;a)`$, is the same as requiring

$$\underset{i=1}{\overset{m}{}}\left(x_ia\right)=0,$$

and any deviations in one will immediately lead to deviations in the other. Hence, they must be proportional to one another. Choosing the coefficient of proportionality, as the second derivative of some appropriate scalar function, $`V`$, gives

$$\frac{}{a}\psi (x_i;a)=V^{\prime \prime }(a)(x_ia),$$
(2)
where the prime stands for differentiation with respect to the argument. The scalar potential, $`V(a)`$, must be independent of the $`x_i`$ because the left-hand side is only a function of $`x_i`$ and a similar equation for $`x_j`$ would lead to a contradiction. We assume that the potential is such that $`V(a^0)=0`$. Consequently, (2) can be rewritten as

$$\frac{}{a}\psi (x_i;a)=\frac{}{a}\left\{V^{}(a)(x_ia)+V(a)\right\}.$$

Integrating from $`a_0`$ to $`a`$ gives

$$\psi (x_i;a)=\psi (x_i;a^0)+\lambda (a)(x_ia)+V(a),$$

where $`\lambda (a)=V^{}(a)`$. In the usual case that the log-likelihood function is logarithmic, we get an exponential family of distributions

$$\mathrm{log}p(x_i;a)=\mathrm{log}q(x_i)+\lambda (a)(x_ia)+V(a)$$
(3)
Averaging both sides with respect to the probability distribution $`P`$ gives

$$\underset{i=1}{\overset{m}{}}\left\{p(x_i;a)\mathrm{log}p(x_i;a)p(x_i;a)\mathrm{log}q(x_i)\right\}=V(a)$$

In information theory, the first term is the negative of the Shannon entropy, the second term is the *inaccuracy*, and the right hand side is the *error* . On the strength of Shannon’s inequality,

$$\underset{i=1}{\overset{m}{}}p(x_i;a)\mathrm{log}\left(\frac{p(x_i;a)}{q(x_i)}\right)=V(a)0$$
(4)
the inaccuracy cannot be smaller than the Shannon entropy. Shannon’s inequality follows very simply from the arithmetic-geometric mean inequality,

$$\underset{i=1}{\overset{m}{}}x_i^{p(x_i;a)}\underset{i=1}{\overset{m}{}}x_ip(x_i;a)$$

with $`x_i=q(x_i)/p(x_i;a)`$.

When $`Q`$ is the uniform distribution, *i.e.*, $`q(x_i)=1/m`$ $`i`$, Shannon’s inequality, (4), becomes

$$S_0(1/m)S_1(P)=V(a)$$

which we have referred to as the entropy reduction caused by the application of a constraint that produces a finite value of $`a`$ . $`S_0(1/m)=\mathrm{log}m`$ is the maximum entropy, and it is known as the Hartley entropy in information theory. Classically, the entropy is defined to within a constant; only entropy differences are measurable.

However, all that we have said so far does not apply to equilibrium thermodynamics . If we average (3) with respect to the $`Q`$ distribution, instead of the $`P`$ distribution, and use Shannon’s inequality, $`_{i=1}^mq(x_i)\mathrm{log}\left[q(x_i)/p(x_i;a)\right]0`$, we immediately run into a problem because $`V(a)`$ must now be necessarily negative. In statistical mechanics, $`q(x_i)`$ represents the surface of constant energy of a hypersphere of high dimensionality . Because of its high dimensionality, the volume of the hypersphere lies very close to its surface so that $`q(x_i)`$ can be thought of as the volume of phase space occupied by the system. Averages are performed with respect to this non-normalizable prior probability distribution . In order to keep the error $`V(a)`$, which will soon be identified as the thermodynamic entropy, positive, it is necessary to introduce a sign change in (3).

This sign change can be rationalized in the following way. The exponential factor, $`e^{\lambda (a)x_i}`$, will not overpower the rapidly increasing factor of the density of states, $`q(x_i)`$. However, the density of states cannot increase faster than a certain power of the radius, $`x_i`$ , of the phase space volume, which is proportional to $`x_i^m`$ in a hypersphere of $`m`$-dimensions. What is needed is an even more rapidly decreasing exponential factor $`\mathrm{exp}\left(\lambda x_i\right)`$.

According to the Boltzmann-Planck interpretation, $`q(x_i)`$ is not a normalized probability, but, rather, a ‘thermodynamic’ probability, being proportional to the volume of phase space occupied by the system. The (random) entropy $`S(x_i)`$ is defined as the logarithm of the thermodynamic probability

$$S(x_i)=\mathrm{log}q(x_i)$$

The phase average is given by

$$E[X]=\underset{i=1}{\overset{m}{}}x_iq(x_i)/\underset{i=1}{\overset{m}{}}q(x_i).$$

The thermodynamic entropy is the phase space average of $`\mathrm{log}\left[q(x_i)/p(x_i;a)\right]`$ *viz*.,

$$S(a)=\frac{_{i=1}^mq(x_i)\mathrm{log}[q(x_i)/p(x_i;a)]}{_{i=1}^mq(x_i)}$$

and its Legendre transform

$$S(a)\lambda (a)a=\mathrm{log}𝒵(\lambda ),$$

defines the logarithm of the generating function, $`𝒵(\lambda )`$. The inaccuracy now appears as the difference between the thermodynamic entropy and the average of the random entropies

$$\underset{i=1}{\overset{m}{}}q(x_i)\mathrm{log}p(x;a)=S(a)\underset{i=1}{\overset{m}{}}q(x_i)S(x_i)/\underset{i=1}{\overset{m}{}}q(x_i)0.$$
(5)
The inequality follows from the facts that $`S`$ increases in the wide sense and is concave. The expectation $`a`$ can be taken either with respect to $`P`$ or $`Q`$. The two averages must necessarily coincide for otherwise there would not be a single general thermodynamics, but rather a “microcanonical thermodynamics” and a separate “canonical thermodynamics” . Taken with respect to $`Q`$, (5) is Jensen’s inequality for a concave function, where the $`Q`$ has positive components but are otherwise arbitrary. Taken with respect to the normalized $`P`$, (5) is the Jensen-Petrović inequality , where $`_{i=1}^mp_i(x_i)x_j`$ for each $`j=1,\mathrm{},m`$. The average of $`m`$ variables is likely to be considerably greater than any of its components. Then, if $`S`$ is increasing,

$$S\left(\underset{i=1}{\overset{m}{}}x_ip(x_i)\right)S(x_j)$$

for $`j=1,\mathrm{},m`$. Multiplying by $`q(x_j)`$ and summing gives back (5). This does not mean that $`S(x_i)/x_i`$ should not decrease: A sufficient condition for $`S(_{i=1}^mx_i)_{i=1}^mS(x_i)`$ is that $`S(x_i)/x_i`$ should decrease.

That $`S(x_i)`$ is an increasing function and $`S(x_i)/x_i`$ decreases, *i.e.*, it is anti-star shaped, are the criteria for inequality attenuation . Fluctuations give rise to inaccuracy (5), and, in their absence a function of the average is equal to an average of the function.

Therefore, if the exponential probability distribution, (3), is to coincide with Gauss’ error law, written in terms of the concavity of the entropy,

$$\mathrm{log}p(x_i;a)=S(x_i)S^{}(a)(x_ia)S(a)=\frac{1}{2}S^{\prime \prime }(\stackrel{~}{a})(x_ia)^2$$
(6)
where $`\stackrel{~}{a}`$ lies between $`x_i`$ and $`a`$, then sign changes are needed. When this is done (3) becomes

$$\mathrm{log}p(x_i;a)=\mathrm{log}q(x_i)\lambda (a)(x_ia)V(a).$$
(7)
A comparison of (6) and (7) shows that the entropy, $`S(a)`$, is the potential, $`V(a)`$, that determines the error law . The concavity of the entropy ensures that the exponent will be negative and hence $`p(x_i;a)`$ will be less than unity. The parameter, $`\lambda (a)`$, is still the derivative of the scalar potential, $`V(a)`$, but since this potential now coincides with the thermodynamic entropy, $`S(a)`$, the Lagrange multiplier $`\lambda (a)`$ is now identified as the internal variable in the entropy representation.

Information theoretic entropies, and the entropy reduction of the thermodynamics of extremes , are not amenable to the previous thermodynamic interpretation, where the entropy is defined as the logarithm of the volume of phase space occupied by the system. Since all the volume lies very near of the surface in a thermodynamic system of high dimensionality, the volume of phase space will coincide with the surface area, which is referred to as the structure function . 111In what turned out to be a futile attempt to justify Tsallis’ formalism, Plastino and Plastino  considered a structure function for the energy of the form $`E_i^{m1}`$. Assuming a bounded phase space—for no given reason— whose total energy is $`E_0`$, they identified the Tsallis exponent as $`\alpha =(m2)/(m1)`$, and, at the same time, *defined* the inverse temperature as $`\beta =(m1)/E_0`$. What they failed to realize is that in order to define a temperature $`m`$ must be much greater than $`1`$ so that $`\alpha 1`$. More precisely $`m`$ must be large enough to validate the use of Stirling’s formula . If the conditions under which they claim Tsallis’ statistical mechanics applies, then it cannot be applied to thermodynamic systems for such systems would be far to small to be capable of defining intensive quantities like temperature and pressure. Rather, we consider the $`P`$ and $`Q`$ as two sets of complete probability distributions. For a given probability distribution, $`Q`$, we seek the set of probability $`P`$ which most closely resemble $`Q`$. This is the minimum discrimination statistic of Kullback .

In order to derive the nonadditive entropies of degree-$`\alpha `$, the logarithm is replaced by the well-known elementary limit

$$\mathrm{log}p(x_i;a)\frac{p^{\alpha 1}(x_i;a)1}{\alpha 1}=\psi (x_i;a),$$

and a similar relation for $`\mathrm{log}q(x_i)`$ in the exponential law we get

$$\frac{p^{\alpha 1}(x_i;a)q^{\alpha 1}(x_i)}{\alpha 1}=\lambda (a)(x_ia)+V(a).$$
(8)
Multiplying (8) by $`p(x_i;a)`$, and summing give

$$_\alpha (Q)S_\alpha (P)=\frac{_{i=1}^mp(x_i;a)\left(p^{\alpha 1}(x_i;a)q^{\alpha 1}(x_i)\right)}{\alpha 1}=V(a)0$$

where

$$_\alpha (Q):=\frac{1_{i=1}^mp(x_i;a)q^{\alpha 1}(x_i)}{\alpha 1}$$

has been referred to as the inaccurary , and

$$S_\alpha (P)=\frac{1_{i=1}^mp^\alpha (x_i)}{\alpha 1}$$
(9)
has been referred to as the Tsallis entropy  in the physical literature, but has been well known in information theory since the late 1960’s . We will henceforth suppress the dependence of the probability distribution $`P`$ on the average $`a`$, because Tsallis’ statistical thermodynamics make no pretext at statistical inference. The inaccuracy is a convex function of $`Q`$, for a given $`P`$, provided $`\alpha 2`$. The inaccuracy is defined as the sum of the entropy of degree-$`\alpha `$, (9), and the error, $`V(a)`$. The inaccuracy has the property that

$$\underset{q1/m}{lim}I_\alpha (Q)=S_\alpha (1/m).$$

The negative of the error is what we have called the entropy reduction, $`\mathrm{\Delta }S`$ .

Now the inequality in (8) follows from Hölder’s inequality. Consider the case when all the $`P`$ are rational; then they can be expressed in the form $`p(x_i)=x_i/_{i=1}^mx_i`$, and $`Q`$ is the uniform distribution, $`q(x_i)=1/m`$ $`i`$. Expression (8) then becomes

$$(\alpha 1)^1\left(\frac{_{i=1}^mx_i^\alpha }{\left(_{i=1}^mx_i\right)^\alpha }m^{1\alpha }\right)=V(a)$$

because of Hölder’s inequalities

$$m^1\underset{i=1}{\overset{m}{}}x_i\stackrel{<}{>}\left(\underset{i=1}{\overset{m}{}}x_i^\alpha \right)^{1/\alpha }m^{1/\alpha }\text{for}\stackrel{\alpha >1}{\alpha <1}.$$

Now, Tsallis and collaborators  find that the maximization procedure of the nonadditive entropy (9) with respect to the constraint

$$a_\alpha =E_\alpha (X)=\underset{i=1}{\overset{m}{}}x_ip^\alpha (x_i)/\underset{i=1}{\overset{m}{}}p^\alpha (x_i),$$
(10)
using the escort probabilities , yields the stationary condition

$$p(x_i)=\frac{[1(1\alpha )\lambda _\alpha (x_ia_\alpha )]^{1/(1\alpha )}}{𝒵_\alpha (\lambda _\alpha )}$$
(11)
where

$$\lambda _\alpha =\lambda /\underset{i=1}{\overset{m}{}}p^\alpha (x_i)$$
(12)
and $`\lambda `$ is the Lagrange multiplier for a constraint, (10). The normalization condition of the $`p(x_i)`$ gives the partition function as

$$𝒵_\alpha (\lambda _\alpha )=\underset{i=1}{\overset{m}{}}[1(1\alpha )\lambda _\alpha (x_ia_\alpha )]^{1/(1\alpha )}.$$
(13)
At best, (11) can be considered as an implicit relation for the probabilities since (12) contains the probabilities explicitly through (12).

In order to reduce Gauss’ principle (8) to something that even vaguely looks like the ‘optimal’ probabilities (11) that maximize the Tsallis entropy, (9), it is necessary to:

1. assume that $`P`$ is an incomplete distribution,
2. set $`q(x_i)=1`$ $`i`$, and
3. set $`V(a)=0`$.

We then obtain

$$\frac{p(x_i)}{_{i=1}^mp(x_i)}=\frac{[1+(\alpha 1)\lambda _\alpha (x_ia_\alpha )]^{1/(\alpha 1)}}{𝒵_\alpha (\lambda _\alpha )},$$

where the partition function is given by (13), and we used the escort probabilities (10) to define the parameter $`a`$, rather than the weighted average (1).

Rather, if we take (7), and introduce the approximation

$$\frac{p^{1\alpha }(x_i)1}{1\alpha }\mathrm{log}p(x_i),$$

and a similar expression for $`q(x_i)`$ we get

$$\frac{p^{1\alpha }(x_i)q^{1\alpha }(x_i)}{1\alpha }=\lambda _\alpha (x_ia_\alpha )V(a_\alpha ).$$
(14)
Setting $`q(x_i)=1`$ and $`V(a_\alpha )=0`$, and requiring the probability distribution $`P`$ to be normalized result in (11), or, equivalently,

$$p^{1\alpha }(x_i)=\frac{1(1\alpha )\lambda _\alpha (x_ia_\alpha )}{𝒵_\alpha ^{1\alpha }(\lambda _\alpha )}.$$

Multiplying both sides by $`p^\alpha (x_i)`$ and summing give

$$𝒵_\alpha ^{1\alpha }(\lambda _\alpha )=\underset{i=1}{\overset{m}{}}p^\alpha (x_i),$$

provided $`\lambda _\alpha `$ is given by the escort average, (10). Rather, if raise both sides of (11) to the power $`\alpha `$, sum, and rearrange, we get

$$𝒵_\alpha ^\alpha (\lambda _\alpha )\underset{i=1}{\overset{m}{}}p^\alpha (x_i)=\underset{i=1}{\overset{m}{}}[1(1\alpha )\lambda _\alpha (x_ia_\alpha )]^{\alpha /(1\alpha )}𝒵_\alpha (\lambda _\alpha ).$$

The only difference between the two forms of averaging is that in the first case use has been made of the escort probability average, (10). Since the two results do not coincide, we conclude that there is something amiss with the escort probability average, (10).

Moreover, if we take the $`\alpha 1`$ limit in (14) we obtain

$$𝒵_1(\lambda _1)=\underset{i=1}{\overset{m}{}}e^{\lambda _1(x_ia)},$$

which is not the correct expression for the partition function even in the unphysical case of a density of states equal to unity.

Finally, multiplying both sides of (14) by $`p^\alpha (x_i)`$, and summing, result in

$$\frac{1_{i=1}^mq(x_i)[p(x_i)/q(x_i)]^\alpha }{1\alpha }=\underset{i=1}{\overset{m}{}}p^\alpha (x_i)V(a_\alpha ).$$

If we now set $`q(x_i)1`$ $`i`$, we come out with

$$S_\alpha (P)=\underset{i=1}{\overset{m}{}}p^\alpha (x_i)V(a_\alpha ).$$
(15)
This shows the correspondence between Shannon entropy and the potential $`V(a)`$ in the $`\alpha 1`$ limit, that was alluded to above in the thermodynamic formulation which takes into account a non-normalized prior probability distribution. However, we have set the prior probability distribution equal to unity, as in the maximum entropy method, and, furthermore, in order to derive the probability distribution (11) from Gauss’ law we had to assume that $`V`$ is identically zero. Relation (15) would, consequently, require the vanishing of the nonadditive entropy, (9).

Based on the foregoing results, we can only conclude that Tsallis’ ‘statistical thermodynamic’ formulation of the nonadditive entropy of degree-$`\alpha `$ is neither correct nor self-consistent.


text:
# Bars and the connection between dark and visible matter

## 1. Introduction

It is impossible to observe dark matter directly, but its existence and a number of its properties can be deduced from its effects on other, visible, galactic components. Thus, properties of the disc and of its substructure can, if correctly interpreted, give us clues on the properties of dark matter haloes. Here I will discuss the connection between bars and dark matter, and what information the former can give us on the latter.

## 2. Angular momentum exchange. Analytical results

The growth of bars in isolated disc galaxies is governed by the exchange of angular momentum between different parts of the galaxy. To understand this better it is best to start with the pioneering paper of Lynden-Bell & Kalnajs (1972; LBK). Using linear analytic theory, these authors showed that it is mainly material at resonance that gains or loses angular momentum. Material at the inner Lindblad resonance (ILR) will lose angular momentum, while material at corotation (CR) and the outer Lindblad resonance (OLR) will absorb it. Since the spiral/bar within CR is a negative angular momentum perturbation, feeding it with angular momentum will damp it, while taking angular momentum from it will excite it. Following in their footsteps, Athanassoula (2003; A03) added a halo (or, more generally, a spheroidal component) and applied the results to the case of bars. Provided the halo distribution function is a function of the energy only, halo material at all resonances will gain angular momentum. This result should be generalisable to other distribution functions, provided energy is the main functional dependence and an appropriate perturbation expansion can be used. Both for the disc and for the halo, there is more angular momentum lost/gained at a given resonance if the density is higher there and if the resonant material is colder and thus more responsive (A03).

So, if a disc galaxy has no halo, or if the latter cannot participate in the angular momentum exchange, the inner disc will emit angular momentum, which will be absorbed by the outer disc. On the other hand, if a halo is present and responsive, then it also will absorb angular momentum. So more angular momentum can be extracted from the inner parts in the presence of a responsive dark matter halo and the bar will grow stronger than in its absence (A03). This explains the ‘paradox’ that bars in halo dominated disc galaxies may grow stronger than in disc dominated cases (Athanassoula & Misiriotis 2002; AM02).

Tremaine & Weinberg (1984) and Weinberg (1985) used nonlinear theory to follow the effect of angular momentum exchange on the slowdown of the bar. They find that, as the bar loses angular momentum, it will slow down, as expected. These works, put together with the analytical studies discussed above, lead to the firm prediction that bars should become stronger and rotate slower in the presence of massive and responsive haloes.

## 3. Angular momentum exchange. Simulation results

In Athanassoula (2002; A02) and A03 I tested the above analytical results with the help of $`N`$-body simulations and I will follow the same path here. The first point to test is that there are indeed stars at (near-) resonance in the halo and that they absorb angular momentum, as predicted by the analytical results. This is shown in Fig. 1, which has been obtained as described in A02 and A03. The upper panels plot the mass per unit frequency ratio, $`M_R`$, as a function of that frequency ratio, namely $`(\mathrm{\Omega }\mathrm{\Omega }_p)/\kappa `$. Here $`\mathrm{\Omega }`$ is the angular frequency, $`\kappa `$ is the epicyclic frequency and $`\mathrm{\Omega }_p`$ is the bar pattern speed. These panels show that the distribution is not at all uniform, but has strong peaks at the resonances. This is true both for the disc, as expected, but also for the halo, as shown initially in A02. The lower panels show the change of angular momentum with time (A03), again as a function of the frequency ratio. Note that the theoretical predictions of LBK, as well as those of A03, are well confirmed by the above results, since these show that a large fraction of both the disc and the halo particles are at resonance, and that they emit/absorb angular momentum, as predicted.

I next check how important the effect of the resonances is. This can be seen in Fig. 2, where I compare the face-on view of the results of three simulations. Initially, the leftmost is disc dominated in the inner parts, or, in the notation of AM02 and A03, it is of MD-type. On the contrary, the simulation shown in the middle panel has initially a large relative halo mass in the inner parts, i.e. it is of MH-type. Finally, the simulation at the right is initially exactly the same as the previous one, except that the halo is rigid (RH). In this case, the halo is only present as a rigid potential, so it can not participate in the angular momentum exchange. The difference between the middle and right panels is indeed stunning! Instead of the strong bar of the MH case, the RH case displays a short oval, confined to the inner parts. Furthermore, comparing the left and middle panels one can clearly see that the strongest bar has grown in the most halo dominated environment. These comparisons, and many others in AM02, A02 and A03, argue strongly for the role of resonances in bar evolution.

Theory predicts that both the amount of mass at resonance and its velocity dispersion should influence the amount of angular momentum exchanged, and therefore the bar strength and the pattern speed decrease. A03, with the help of appropriate sequences of simulations, showed this to be indeed true in the case of $`N`$-body simulations. Thus, the increase of the bar strength and the decrease of its pattern speed are jointly set by three factors : the relative amount of halo mass at (near-) resonance, how hot the disc is at the resonances, and how hot the halo is at the resonances. Any one on its own is not sufficient to determine the outcome, since all three can limit the amount of angular momentum exchanged.

Since the amount of angular momentum exchanged is a determining factor, one expects to find correlations between this quantity and the bar evolution. Unfortunately, the angular momentum exchanged is a difficult quantity to measure, so in Fig. 3, as in A03, I use instead the total angular momentum gained by the halo, relative to the total angular momentum initially in the system. This, however, is a good approximation of the total angular momentum exchanged only in cases when the outer disc absorbs only little angular momentum. Even so, Fig. 3 shows, for a very large number of simulations, a clear correlation between the relative amount of angular momentum absorbed by the halo and the bar strength, as well as a trend between it and the bar pattern speed. Furthermore, if I limit myself to simulations where the halo is the main absorber (lower panels), then the latter trend becomes a very strong correlation (middle lower panel). These results are basically the same as those shown in A03, except that I have added here a few more simulations, which have been run since that paper was written.

How much are the above correlations dependent on the particular type of initial conditions used? It is reasonable to assume that initial conditions having different density and velocity profiles for the disc, bulge and halo component will also give such correlations, since these correlations are the reflection of the physics underlying bar evolution. However, it could well be that the regression lines (or loci defined by the trend) are somewhat differently located (shifted) on the corresponding planes. That can be properly tested only by repeating this type of work with different initial models; a rather daunting task seen the very large number of simulations necessary. Yet a few clues exist already. For example, from the middle lower panel one can see that the regression line is somewhat shifted to the right for simulations with a more concentrated spheroidal. Similarly, the right bottom panel shows that the regression line is somewhat lowered for simulations with initially hotter discs. The differences, however, are rather small, and one can conclude that more angular momentum exchange will lead to stronger bars that rotate slower.

Berentzen et al. (2003) present a study where the bar is driven by a companion in a disc in which a previous bar had been destroyed by gas inflow. They find good agreement of the results of their simulations with those presented here on the ($`S_B`$, $`\mathrm{\Omega }_p`$) plane, even though the initial models and particularly the problem at hand are very different.

How does the amount of angular momentum exchanged influence the morphology of the bar? Fig. 4 shows the results of three simulations. In the one shown in the left panel a large amount of angular momentum was exchanged during the evolution. It has a strong bar, as is the case for all such simulations. The middle and right panels show examples from simulations where little angular momentum has been exchanged. In the middle one, the disc is initially hot and the halo is relatively cool and concentrated. Thus the outer disc can not contribute much to the angular momentum exchange and it is the halo that absorbs most of the angular momentum. The result is an oval, which extends to large radii, in order to maximize the emitting region. In the last example (right panel) the halo is very hot, so it can not absorb much angular momentum, while the disc is initially relatively cool. Thus the bar has to stay short, in order for its CR and OLR to be in regions of relatively high density, so as to contribute an efficient angular momentum sink. To summarize, one can say that cases where a lot of angular momentum has been exchanged always display strong bars. On the other hand, in cases where only little angular momentum has been exchanged, the bar can either be short, or longer but more like an oval, depending on whether it is the disc or the halo that is the main absorber.

## 4. Why is there such a wide variety of bar strengths in observed galaxies, and why are some galaxies non-barred ?

Observed bars come in a wide variety of strengths, ranging from very strong bars, like NGC 4608 (Gadotti & de Souza 2003) or NGC 1365, to short bars, like our own Galaxy, or ovals, like in NGC 1566. Furthermore, a few discs have no bar at all (Grösbol, Pompei & Patsis 2002). How does that compare with the $`N`$-body results?

Haloes were initially thought to stabilise discs (e.g. Ostriker & Peebles 1973). More recent work (AM02, A02, A03) shows that this statement is not necessarily correct, since haloes can stimulate bar growth by taking angular momentum from the inner disc. Nevertheless, bars forming in disc galaxies with a substantial halo take longer to form than bars forming in galaxies with less halo, so that, in a sense, a halo can be considered as having a stabilising influence, although the bars that grow in a more massive and responsive halo may finally become stronger. A03 showed that this can happen for a wide range of halo to disc mass ratios, provided of course the halo is responsive. This, however, will not necessarily extend to very low relative disc masses (A03). Simulations to test this limit would be very CPU intensive, since they would necessitate a very large number of particles, and also because the bar growth would be exceedingly slow. Moreover, if the bar grows in time scales much longer than the Hubble time, the problem is of little astronomical interest and the corresponding galaxies can be considered bar-stable.

Haloes can have a further stabilising influence either if their mass distribution is such that not much material is at resonance with the bar, or, more likely, if the resonant halo material is very hot and thus can not absorb much angular momentum. Unfortunately, not much is known about the velocity distribution within the halo component in real galaxies.

Similarly, a bar forming in a hot disc will take longer to grow than in a cold one (Athanassoula & Sellwood 1986, A03). In very hot discs, when the bar eventually forms, it has the form of an oval (see Fig. 4).

Thus, with the help of $`N`$-body simulations, we can account for the large variety of bar strengths observed in real galaxies. As discussed also in the previous section, this will be determined by the total amount of angular momentum exchanged and also by the specific amount by which each of the partners (inner disc, outer disc, halo, bulge) enters in the exchange. In this picture, galaxies with no bars should have a very low relative disc mass and/or a very hot disc and/or a rather unresponsive halo, so that the bar takes very long to grow and so that the halo can not help its secular evolution.

This picture is not complete. There is still the question of black holes, and other ‘extreme’ central concentrations, which should have a stabilising influence. The interplay between this and the responsiveness of the disc and halo will be discussed elsewhere.

## 5. A bar in the halo

In cases with a sufficiently strong bar in the disc component, the halo does not stay axisymmetric, but shows also a bar-like or oval deformation. An example can be seen in Fig. 2 of Holley-Bockelmann, Weinberg & Katz (2003; HBWK). I found such structures also in my own simulations and I call them, for simplicity, halo bars. Preliminary results show that, in cases with a strong disc bar, the halo bar is triaxial, but prolate like, with axial ratios of the order of 0.7 or 0.8 in the inner parts, and becomes more spherical as the radius increases. The length of the halo bar increases with time, but always stays considerably shorter than that of the disc bar. It is roughly aligned with the disc bar at all times (at least within the measuring errors), i.e. it turns with roughly the same pattern speed. This means that it is a slow bar. The bisymmetric component in the halo extends well beyond the end of the halo bar, and there it trails behind the disc bar, much as seen in Fig. 2 of HBWK. The properties of these halo bars will be discussed in detail elsewhere (Athanassoula, in prep.).

## 6. Do bars flatten or steepen central halo cusps?

Current cosmological CDM simulations predict that dark matter haloes are cuspy (e.g. White, this volume), although there is still no agreement reached about the value of the inner slope. On the other hand, observations have shown that, at least in a number of cases where the data are of sufficiently good quality, haloes have cores (e.g. Bosma, or de Blok, this volume). If haloes were indeed formed with a cusp and now have a core, then some mechanism during galaxy formation and/or evolution should be responsible for this change. Several have been so far proposed, of which one involves a bar. Indeed, Hernquist & Weinberg (1992) and Weinberg & Katz (2002), based on analytical calculations and on $`N`$-body simulations of a cusped halo containing a rigid bar and no disc, proposed that the bar could flatten the cusp by giving angular momentum to the halo. This, however, was not generally accepted, since the fully self-consistent simulations of Sellwood (2003; S03) and Valenzuela & Klypin (2003; VK) showed the contrary, i.e. that during the evolution the halo radial profile showed a small, yet clear, steepening. These authors attribute the flattening found by the previous simulations to the fact that they do not include a live disc, and thus neglect the slow gradual change of the resonant positions due to the slowdown of the bar, and the disc responsiveness. This explanation did not, in turn, satisfy HBWK, whose fully self-consistent simulations showed a clear flattening of the halo cusp (HBWK) and thus argued that self-consistency, or the lack of it, could not be the factor determining why some of the previous simulations showed a flattening and others a steepening. HBWK instead argue that the cusp steepening found by S03 and VK is spurious, and due to the fact that their simulations are noisy and thus do not describe the resonances sufficiently accurately. Numerical noise could indeed, if sufficiently strong, prevent the halo resonances from absorbing the correct amount of angular momentum and thus lead to a wrong evolution. Other numerical shortcomings, however, e.g. leading to a wrong coupling between the planar and vertical resonances, could also influence the results, producing a spurious flattening or steepening.

Faced with the difficult task of reviewing this unsatisfactory state of affairs, I turned to my own simulations (AM02, A02, A03) to see whether they could provide any clues. They are neither grid based, nor use SCF, and thus could provide an independent view. They were not specifically designed to tackle this problem, so the initial halo has a core (albeit sometimes small). They can, nevertheless, lead to a number of insightful conclusions, which I will discuss briefly below. In particular, my study of the resonances allows me to single out and to follow the resonant stars. I could thus make sure that they cover the phase space adequately. Furthermore, in A02 I showed that there is a large trapping of the particles at resonances, which shows that they are not unduly knocked off their trajectories by noise during the evolution. Thus, my simulations fulfill all the necessary conditions set by HBWK to perform the test at hand.

In order to check whether the density profile in the inner parts flattens or steepens, I simply measured the amount of mass in concentric spherical shells111Although very straightforward, this method is not the most appropriate for the problem at hand, because, as discussed in the previous section, the halo isodensities are triaxial and not spherical. Results based on other methods of calculating the density will be discussed elsewhere.. Since the existence of a bar is necessary for the mechanism to work, I confined myself to times after the bar had grown. I find that, in the vast majority of the cases I checked, there is a steepening of the halo radial density profile, albeit not large. In a couple of cases, however, I noted a very small flattening. Since this is even smaller in amplitude than the small steepening found in the remaining cases, I was ready to discard it as insignificant, until I noticed that it occurred in the more centrally concentrated cases. This prompted me to investigate the problem further.

The analytical predictions are clear and follow directly from what was discussed in the previous sections : Halo material at (near-) resonance should absorb angular momentum and move to larger cylindrical radii222It is important in this problem to distinguish between cylindrical and spherical radius.. For the case of the ILR, such material is located in the inner parts of the halo and this should lead to a flattening of the cusp. I have already shown in the previous sections that my simulations confirm the theoretical predictions about the existence of the resonances and the angular momentum exchanged. Showing the increase of the cylindrical radius of particles on orbits at, or around, the ILR is no an easy task, since they are trapped around periodic orbits of the $`x_1`$ tree333The $`x_1`$ tree is the 3D orbital analog of the 2D $`x_1`$ family and provides the backbone of the 3D bars (Skokos, Patsis & Athanassoula 2002)., for which it is not always easy to calculate the frequencies and the time average of the radius. Moreover, there are relatively few particles at ILR in the models I simulated (Fig. 1). On the other hand, it is very easy to see the increase in cylindrical radius in the case of CR, which does not suffer from the above drawbacks.

This, however, is not the complete picture, since there are three effects working against the previous one, and thus leading to a steepening of the halo :

* As stressed by AM02, VK and O’Neill & Dubinski (2003), the radial density profile of the disc becomes considerably more centrally concentrated with time. In so doing, the disc material pulls the halo material also inwards, thus causing it to contract. This leads to an increase of the halo radial density profile in the inner parts.
* As the halo gains angular momentum it will get flattened towards the disc equatorial plane and this will result in a decreasing of the spherical radius of the individual particles. As seen in the previous section, this happens mainly in the inner parts of the halo and will, therefore, lead to a steepening of the cusp.
* During the evolution the bar becomes stronger and this influences the shape of the periodic orbits of the $`x_1`$ tree and of the regular orbits trapped around them. Namely, their axial ratio $`a/b`$ increases and they approach nearer to the center, so that sometimes the trapped orbits can actually cover the central area. This means that, although their average cylindrical radius may increase, the central-most area may have an increased density.

There are thus competing effects : On the one hand the resonant particles can move to larger cylindrical radii and thus tend to flatten the cusp. On the other hand, there are other effects, also linked to the angular momentum redistribution, which tend to diminish the spherical radii of the particles. It is the outcome of the competition between these effects that will decide whether the inner density profile will become steeper or shallower during the evolution. It is thus not necessarily surprising that HBWK find a flattening of the cusp, while S03 and VK find a steepening. The final result will depend on the distribution function of the halo, of which hardly anything is known, as well as that of the disc. In as far as the $`N`$-body simulations are concerned, it could also, unfortunately, depend on numerical effects, which could artificially modify the effect of the resonances.

The application of this mechanism to real galaxies faces, to my mind, two further serious problems. One is that haloes have substructure. Although this will not inhibit resonance driven evolution (Weinberg 2001), it could still perturb the effect of the resonances. The second, perhaps even more serious, problem is the presence of gas. All simulations that have so far studied this mechanism are purely stellar and have no gas. A gaseous component could increase further the central concentration of the disc material (e.g. Athanassoula 1992, Heller & Shlosman 1994) and thus make its inwards pull on the halo material yet stronger. It would then become yet more difficult for the halo (near-) resonant stars to overcome the inwards pull and achieve a flattening of the cusp.

## References

Athanassoula, E. 1992, MNRAS, 259, 345

Athanassoula, E. 2002, ApJ, 569, L83 (A02)

Athanassoula, E. 2003, MNRAS, 341, 1179 (A03)

Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35 (AM02)

Athanassoula, E., & Sellwood, J. A. 1986, MNRAS, 221, 213

Berentzen, I., Athanassoula, E., Heller, C. H., & Fricke, K. J. 2003, astro-ph/0309664

Gadotti, D. A. & de Souza, R. E. 2003, ApJL, 583, L75

Grösbol, P., Pompei, E. & Patsis, P. A. 2002, in “Galactic Discs : Kinematics, Dynamics and Perturbations”, eds. E. Athanassoula, A. Bosma and R. Mujica, PASP conference series, 275, 305

Heller, C. H., & Shlosman, I. 1994, ApJ, 424, 84

Hernquist, L., & Weinberg, M. D. 1992, ApJ, 400, 80

Holley-Bockelmann, K., Weinberg, M. D., & Katz N. 2003, astro-ph/0306374 (HBWK)

Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1, 1972 (LBK)

O’Neill, J. K., & Dubinski, J. 2003, MNRAS, astro-ph/0305169

Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467

Sellwood, J. A. 2003, ApJ, 567, 638 (S03)

Skokos, Ch., Patsis, P., & Athanassoula, E. 2002, MNRAS, 333, 847

Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729

Valenzuela, O., & Klypin, A. 2003, MNRAS, 345, 406 (VK)

Weinberg, M. 1985, MNRAS, 213, 451

Weinberg, M. 2001, MNRAS, 328, 321

Weinberg, M., & Katz, N. 2002, ApJ, 580, 627


text:
# Directed diffraction without negative refraction

## Figure Captions

(a) The band structure of a square lattice of air holes embedded in the dielectric medium $`ϵ=12`$. The lattice constant is $`a`$ and the radius of the cylinders is $`0.35a`$. $`\mathrm{\Gamma }M`$ and $`\mathrm{\Gamma }X`$ denote the  and  directions respectively. (b) The magnetic field $`H_z`$ of a point source and its image across a photonic crystal slab called superlens in Luo et al.PRB . The slab measures 11.31x31.11, and the frequency is 0.192. (c) The conceptual layout of a superlensing phenomenon, adapted from VeselagoVes  and PendryPendry . (d) The illustration of negative refraction of a plane wave across a flat slab, by analogy with the right panel of Fig. 2 in Ref. PRB .

The intensity image of the magnetic fields across a flat photonic crystal slab made of air holes in dielectric $`ϵ=12`$. The principle directions of the crystal are drown in the figure, as  and . A plane wave is incident to the slab with various incident angles in degrees with respect to the  direction: (a) 0; (b) 22.5; (c) 45; and (d) 67.5. The slab measures as 49.5x92 in terms of lattice constant. For clarity, the air holes are not shown in the figure.

The intensity image of the magnetic fields across a flat photonic crystal slab. The principle directions of the crystal are drown in the figure, as  and . The  direction is tilted towards left by an angle of 22.5 degree with reference to the normal direction of the slab surface. A plane wave is incident to the slab with various incident angles in degrees with respect to the normal direction of the slab surface: (a) 0; (b) 22.5; (c) 45; and (d) 67.5. The slab measures as 49.5x92 in terms of lattice constant.

The intensity image of the magnetic fields across two flat photonic crystal slabs. The principle directions of the crystal are drown in the figure, as  and . The two slabs measure in terms of lattice constant as (a) 70x92;(b) 49.5x120. A plane wave is incident to the slab with the incident angle of 22.5 degree with respect to the normal direction of the slab surface.


text:
# Formation of star clusters in the LMC and SMC. I. Preliminary results on cluster formation from colliding gas clouds

## 1 Introduction

It well established that several physical properties of the globular clusters and populous young blue clusters in the Large Magellanic Cloud (LMC) are differ markedly from those of clusters in the Galaxy (e.g., van den Bergh 2000a). These properties include the more flattened shapes of the LMC clusters (e.g., Geisler & Hodge 1980; van den Bergh & Morbey 1984), the disky distribution of its globular cluster system (e.g., Schommer et al. 1992), a larger fraction of apparently binary clusters or physical cluster pairs in the LMC (Bhatia & Hatzidimitriou 1988; Bhatia et al. 1991; Dieball & Grebel 1998), a possible “age/metallicity gap” (e.g., Da Costa 1991; Olszewski et al. 1991; Geisler et al. 1997; Sarajedini 1998; Rich et al. 2001), larger sizes at a given galactocentric distance (van den Bergh 2000b), and the presence of a significant number of massive young to intermediate age clusters in the LMC (e.g., van den Bergh 1981, 2000a).

The higher fraction of binary clusters in the LMC, in particular, has attracted much attention from theoretical and numerical works. Fujimoto & Kumai (1997) proposed that oblique cloud-cloud collisions in an interaction between the LMC and the Small Magellanic Cloud (hereafter, SMC) can result in the formation of binary star clusters revolving around each other. Leon et al. (1999) proposed a different scenario in which the tidal capture of clusters in a group (where tidal encounters are expected to be more common) could be associated with the formation of the LMC binary clusters. de Oliveira et al. (2000) suggested that merging of binary clusters can be responsible for the observed flattened shapes of LMC clusters.

Kumai et al. (1993) and Fujimoto & Kumai (1997) pointed out that if interstellar gas clouds are in large-scale disorganized motions with velocities of more than $`50100`$ km s-1 in the interacting LMC/SMC system, then they may collide with one another to form compact star clusters due to strong shock compression. This idea is supported by the coincidence between the observationally inferred two “burst” epochs ($``$ 100 Myr and $`12`$ Gyr ago of cluster formation after the initial LMC collapse phase) (e.g., Girardi et al. 1995) and the theoretically predicted epochs of closest encounter between the LMC and the SMC (Murai & Fujimoto 1980; Gardiner et al. 1994; Gardiner & Noguchi 1996). However, due to the lack of extensive numerical studies of this scenario, these authors did not address (1) whether the high-speed, oblique cloud-cloud collisions, which are central to the scenario of Kumai et al. (1993), are present in the LMC/SMC interaction and (2) how star formation efficiency may increase in the colliding clouds such that compact and bound star clusters rather than unbound field stars will be formed.

In this paper, by using numerical simulations, we demonstrate that the star formation efficiency of colliding gas clouds in interacting galaxies can significantly increase, resulting in the formation of compact stellar systems. This numerical investigation is two-fold: We first derive the most probable impact parameter and the relative velocity of two colliding clouds in a large-scale dynamical simulation of the interacting LMC/SMC. We then investigate the hydrodynamical evolution and star formation processes of colliding clouds, based on the most probable parameter values for cloud-cloud collision derived in our first set of simulations. Here we describe the formation of star clusters in colliding clouds in a general way, rather than attempt to explain precisely the observed physical properties of star clusters for the LMC/SMC system in a self-consistent manner. Since this paper is the first step toward the better understanding of star cluster formation in the LMC-SMC system, the numerical models are somewhat idealized and lack realism in some areas. In future papers, we will describe the formation processes, structure and kinematics, and chemical properties of star clusters formed from cloud-cloud collisions, using a more sophisticated simulation. The origin of disky distribution of LMC old star clusters will be also discussed in our future papers in terms of recent cosmological simulations of globular clusters (e.g., Kravtsov & Gnedin 2003).

## 2 Models

We first describe the disk galaxy models in which the dynamical evolution of gas clouds in interacting LMC/SMC is investigated. We then describe the hydrodynamical models for the evolution of colliding gas clouds with star formation in the LMC/SMC system.

### 2.1 Tidal interaction in the LMC/SMC system

We model LMC and SMC as bulge-less, gas-rich disks with an initial gas mass fraction of 0.1. They Magellanic Clouds are modeled in a fully self-consistent way by using the Fall-Efstathiou (1980) model, with the exponential density profile for the disks and halo dark matter-to-disk mass fraction equal to 4. The total galactic mass and the disk size for the LMC (SMC) are assumed to be 2.0 $`\times `$ $`10^{10}`$ $`\mathrm{M}_{}`$ (2.0 $`\times `$ $`10^9`$ $`\mathrm{M}_{}`$) and 7.5 kpc (2.4 kpc), respectively. In order to investigate the nature of cloud-cloud collisions in the interacting LMC/SMC, we adopt the “sticky particle method” (e.g., Hausman & Roberts 1984) in which the interstellar medium is described as an ensemble of discrete gas clouds. The size ($`r_{\mathrm{cl}}`$) of an individual cloud with a given mass ($`M_{\mathrm{cl}}`$) is chosen such that the cloud satisfies the observed mass-size relation of gas clouds (Larson 1981). All calculations related to self-gravitating gas clouds and stellar components were carried out on GRAPE systems (Sugimoto et al. 1990) and total particle number in each simulation is 20000 for dark matter and 22000 for disk components.

We focus on the past 1 Gyr evolution of the LMC/SMC (in particular, at the latest SMC’s pericenter passage), during which time populous young star clusters are known to have formed (e.g., van den Bergh 2000a). Since the tidal effect on the LMC/SMC system due to the Galaxy is small compared to that from the interaction between the Magellanic Clouds themselves at SMC’s pericenter passage (Gardiner & Noguchi 1996), we do not explicitly include the gravitational effect of the Galactic dark matter halo. Guided by the earlier numerical results of tidal interaction between the Galaxy, the LMC, and the SMC (Gardiner et al. 1994; Gardiner & Noguchi 1996), we choose a plausible set of orbital parameters for the LMC/SMC. The apocenter radius of the interaction (for the last 1 Gyr) is set to be 30 kpc. The pericenter of the orbit (represented by $`R_\mathrm{p}`$), inclination of the LMC’s disk with respect to the orbital plane ($`\theta _{\mathrm{LMC}}`$), and that of SMC ($`\theta _{\mathrm{SMC}}`$) are assumed to be free parameters. Although we investigated a number of models with different parameter values, we present the results of the model with $`R_\mathrm{p}`$ = 3.75 kpc, $`\theta _{\mathrm{LMC}}`$ = $`15^{}`$, and $`\theta _{\mathrm{SMC}}`$ = $`45^{}`$. We choose this set of parameters since they exhibit behavior characteristic of cloud-cloud collisions (e.g., distribution of relative velocity of clouds) in the present study.

The most important parameter of this model, involving interacting galaxies with a mass ratio of 0.1, is the pericenter of the SMC. If the pericenter distance is large (e.g. 15 kpc, which is twice the LMC’s disk size), the frequency of cloud-cloud collisions is not enhanced significantly in the simulation, and the star formation rate is not significantly increased in the simulated LMC disk. Accordingly, the parameter of the pericenter should be carefully chosen. We base the parameter values of the LMC/SMC orbital properties on the early numerical simulations by Gardiner & Noguchi (1996), which are not only consistent with the observed location and radial velocity of the LMC/SMC system but also successful in reproducing the observed physical properties of the Magellanic stream. The orbital parameters in the present study (and thus in Gardiner & Noguchi 1996) are broadly consistent with $`HIPPARCOS`$ data by Kroupa & Bastian (1997), as shown by Yoshizawa & Noguchi 2003. Therefore, the adopted parameter values can be regarded as reasonable, though the exact orbital properties of the LMC/SMC system have not been observationally determined (thus the orbital parameters should be still free parameters).

It is possible that introducing a mass spectrum on clouds instead of a single mass for all clouds affects the results shown in this paper. The total number of cloud-cloud collisions ($`N_{\mathrm{cl}}`$) between the clouds with masses of $`M_{\mathrm{cl}}`$ during the time interval $`dt`$ can be written as, $`N_{\mathrm{cl}}=\pi \times dt\times r_{\mathrm{cl}}^{}{}_{}{}^{2}\times V_{\mathrm{cl}}\times \mathrm{\Phi }_{\mathrm{cl}}`$, where $`r_{\mathrm{cl}}`$, $`V_{\mathrm{cl}}`$, and $`\mathrm{\Phi }_{\mathrm{cl}}`$ are the cloud radius, typical cloud velocity, and number density of gas clouds. By assuming that the number density for a given volume is proportional to the cloud number function observed in the Galaxy (e.g., Harris & Pudritz 1994) and adopting the Larson’s mass-size relation (1981), we can derive the $`M_{\mathrm{cl}}`$ dependence of $`N_{\mathrm{cl}}`$. Since $`V_{\mathrm{cl}}`$ is highly likely to be independent on $`M_{\mathrm{cl}}`$, the observed relations $`\mathrm{\Phi }_{\mathrm{cl}}M_{\mathrm{cl}}^{}{}_{}{}^{1.63}`$ and $`M_{\mathrm{cl}}r_{\mathrm{cl}}^{}{}_{}{}^{2}`$ imply that $`N_{\mathrm{cl}}M_{\mathrm{cl}}^{}{}_{}{}^{0.63}`$. This derived relation suggests that (1) cloud-cloud collision rates are larger for the clouds with smaller masses and (2) larger clouds can more frequently collide with smaller clouds. Therefore a spectrum of cloud masses is expected to lead to a larger number of low mass clusters.

### 2.2 Star formation in colliding gas clouds

Next we investigate the hydrodynamical evolution of two colliding clouds by using a TREESPH code with star formation (Bekki 1997). The initial cloud mass ($`M_{\mathrm{cl}}`$) and size ($`r_{\mathrm{cl}}`$) are set to be $`10^6`$ $`M_{}`$ and 97 pc, respectively, which are consistent with the observed mass-size relation by Larson (1981), and therefore with the large-scale simulations described in §2. A gas cloud is assumed to have an isothermal radial density profile with $`\rho (r)1/(r+a)^2`$, where $`a`$ is the core radius of the cloud and set to be $`0.2r_{\mathrm{cl}}`$. An isothermal equation of state with a sound speed of $`c_\mathrm{s}`$ is used for the gas, and $`c_\mathrm{s}`$ is set to 4 km s-1 for models with $`M_{\mathrm{cl}}`$ = $`10^6`$ $`M_{}`$. We choose this value of $`c_\mathrm{s}`$ guided by the virial theorem and Larson’s mass-size relation. A gas particle in a given cloud is converted into a collisionless stellar particle if two conditions are met: First, the local dynamical time scale (corresponding to $`(4\pi G\rho _i)^{0.5}`$, where $`G`$ and $`\rho _i`$ are the gravitational constant and the density of the gas particle, respectively) must be shorter than the sound crossing time (corresponding to $`h_i/c_\mathrm{s}`$, where $`h_i`$ is the smoothing length of the gas) Secondly, the gas flow is converging. This method therefore mimics star formation due to the Jeans instability in gas clouds.

The initial orbital plane of the two colliding clouds with relative velocities of $`V_\mathrm{r}`$ and impact parameter of $`b`$ is set to be the $`x`$-$`y`$ plane. The position and the velocity of each cloud is represented by $`\stackrel{}{x}_i`$ and $`\stackrel{}{v}_i`$ ($`i`$ = 1, 2), respectively. We generally only show our “standard model” in which $`\stackrel{}{x}_1`$ = ($`1.5r_{\mathrm{cl}}`$,$`0.5b`$,0) = $`\stackrel{}{x}_2`$, $`\stackrel{}{v}_1`$ = ($`0.5V_r`$,0,0) = $`\stackrel{}{v}_2`$, $`V_\mathrm{r}`$ = 20 km s-1 (or $`V_\mathrm{r}/c_\mathrm{s}`$ = 5), and $`b`$ = $`0.5r_{\mathrm{cl}}`$ (= 48.5 pc), since this model describes the typical behavior of star cluster formation in colliding clouds in our simulations. Using the most probable values of $`V_\mathrm{r}`$ and $`b`$ derived from our large-scale simulations, we investigate the parameter dependencies of formation processes of star clusters on $`V_\mathrm{r}`$ and $`b`$ for 0 $``$ $`V_\mathrm{r}`$ $``$ 67 km s-1 and 0 $``$ $`b/2r_{\mathrm{cl}}`$ $``$ 1. 20000 SPH particles are used in a simulation.

Because our adopted total particle number is limited, the resolution of the simulation is at most $`10^2`$ $`M_{}`$ in mass and $``$ 1 pc in scale for the models with $`M_{\mathrm{cl}}`$ = $`10^6`$ $`M_{}`$. Therefore, a stellar particle converted from a gas particle does not represent directly an individual star with a mass and size the same as that observed. Most of the stars in our simulations are formed in the very center of a gas cloud, where the Jeans mass of the gas is of order $`10^210^3`$ $`M_{}`$ due to the lower temperature and the higher gas density. Therefore, the stellar particles in our simulations can be regarded as small “sub-clusters” of stars with the mass of $`10^210^3`$ $`M_{}`$ Our future higher-resolution simulations with the total (gaseous) particle number of more than $`10^6`$ will enable us to address not only the physical properties of the sub-clusters, but also the formation of a single cluster from the merging between these sub-clusters.

We choose an initial gas temperature (sound speed for isothermal gas) guided by the virial theorem for a gas cloud with a given mass and size. In other words, the value of $`c_\mathrm{s}`$ is chosen such that an isolated gas cloud (not merging with another cloud) is unable to collapse spontaneously. Accordingly, the isolated gas model cannot form stars in its central regions, and the fragmentation of such gas clouds does not occur. This ensures that if star formation occurs in colliding clouds, it is purely a result of the hydrodynamical evolution of gas driven by cloud-cloud collisions in our numerical study. Thus, the adopted assumption of isothermal equation of state and a higher initial gas temperature (corresponding to the virial temperature of the gas) can help us to better interpret the derived results of our numerical simulations. However, the reader should note that this prescription is not as physically realistic as including heating and cooling sources in the gas.

## 3 Results

### 3.1 The probability of high-speed, oblique cloud-cloud collisions

Figure 1 shows how the frequency of cloud-cloud collisions can be enhanced and, which type of cloud-cloud collisions most frequently occur, during the LMC/SMC interaction. As the SMC passes by the pericenter of the orbit at $`T`$ = 0.61 Gyr (where T represents the time that has elapsed since the two disks begin to interact), the strong tidal force induces non-axisymmetric structures (i.e., bars and spiral arms) in the disks of both the LMC and SMC. The large-scale tidal force randomizes the motion of the clouds during the interaction, and consequently the cloud-cloud collision rate is increased by a factor of $`>`$ 4 after the pericenter passage. The cloud-cloud collisions with $`V_\mathrm{r}`$ $``$ $`2540`$ km s-1 are the most common during the interaction (e.g., $`T`$ = 0.92 Gyr), and we estimate the mean $`V_\mathrm{r}`$ to be 60 km s-1. The impact parameter represented by $`b`$ (Binney & Tremaine 1987) for cloud-cloud collisions can be widely distributed during the tidal interaction, although collisions with larger values ($`b/2r_{\mathrm{cl}}`$ $`>`$ 0.5) represent $``$ 2/3 of the total. These results demonstrate that high-speed ($`V_\mathrm{r}`$ $`>`$ 50 km s-1), oblique collisions between two similar clouds are likely in the LMC/SMC interaction, thus confirming the earlier suggestions by Fujimoto & Kumai (1998).

However, it should be stressed here that colliding gas clouds with relative velocities of more than $`50100`$ km s-1 are not particularly common amongst colliding gas clouds in our simulation. Therefore, the proposed large-scale motions with velocity of more than $`50100`$ km s-1 (Kumai et al. 1993; Fujimoto & Kumai 1997) are less likely in interacting LMC/SMC. The reason for the lower velocities here ($`V_\mathrm{r}`$ of $``$ $`2540`$ km s-1, rather than $`50100`$ km s-1) is that the LMC’s cloud system is not strongly disturbed by the SMC tidal field due to the small mass ratio of the SMC to the LMC ($``$ 0.1). Our numerical results shown in Figure 1 suggest that if star clusters are formed from cloud-cloud collisions in interacting LMC/SMC, then the clusters formed from clouds with relative velocities of more than $`50100`$ km s-1 constitute only a minor population amongst the ensemble of young clusters in the LMC/SMC systems. The pros and cons of the original collisional formation model of star clusters in the LMC/SMC system (Kumai et al. 1993) are discussed in detail later.

### 3.2 The formation of star clusters in colliding gas clouds

Figures 2 and 3 illustrate how star clusters are formed in two colliding clouds in our standard model. Strong gas compression and dissipation during the collision leads to an elongated slab-like structure formed at around $`T`$ = 17.1 Myr, where $`T`$ represents the time elapsed since the two clouds began to collide. As the dissipative merging proceeds, the density of gas becomes very high in the shocked regions which are originally the central regions of the two clouds ($`T`$ = 22.8 Myr). Two compact clusters are formed in these high-density gas regions and begin to orbit each other ($`T`$ = 22.8 Myr). This result implies that the orbital angular momentum of the two gas clouds is efficiently converted into that of the binary star clusters during the dissipative cloud-cloud collision. The star formation is akin to an instantaneous “starburst” with a maximum star formation rate of 0.095 $`M_{}`$ yr-1, and 40 % of the gas is converted into stars within 10 Myr.

A stellar particle is assumed to be formed from each gas cloud which is considered to be collapsing due to gravitational instability (i.e., Jeans instability in the present study). The Jeans mass ($`M_\mathrm{J}`$) of gas in the central regions of colliding clouds is estimated to be $`10^3`$ $`M_{}`$, and, as such, each stellar particle can be regarded as representing a small “sub-cluster”. Since all of these sub-clusters are formed in the very center of each colliding cloud, the single massive cluster formed in each colliding cloud may be thought of as consisting of numerous small sub-clusters in the early stages of massive cluster formation. Unfortunately, because of the limited resolution of the present simulation, we cannot investigate the subsequent dynamical evolution of these sub-clusters. However, we may reasonably assume that these numerous sub-clusters will finally merge with one another and consequently erase all of the substructures inside the cluster. In fact, we do not observe any new stellar particles escaping from the parent clouds because they are initially in the deepest potential well (i.e., the center of the clouds). Thus a single star cluster with a very smooth and homogeneous mass distribution would finally form.

The identification of such substructure in our simulations is interesting when one considers recent evidence for the presence subclustering in young cluster environments. Testi et al. (2000) identified three spatially and kinematically distinct subclusters within Serpens, a nearby ($``$300pc) cloud comprising of 500-1500 M$`\mathrm{}`$ molecular gas and 40–80M$`\mathrm{}`$ young stellar objects (Giovannetti et al. 1998). Our simulations support the idea that the characteristic Jeans scale in cloud-cloud collisions gives rise to several subclusters which then proceed to merge in a bottom-up (hierarchical) fashion. The idea that star formation in stellar clusters proceeds preferentially in ”sub-clusters” of enhanced stellar density was suggested by Clarke, Bonnell & Hillenbrand (2000).

The formation process described in our simulations differs in important ways from that proposed originally by Kumai et al. (1993). Kumai et al. (1993) envisaged the formation of star and globular clusters from gravitationally unstable gas with $`M_\mathrm{J}`$ of $`10^510^6`$ $`M_{}`$ in colliding clouds. Our simulations suggest that a single massive cluster is initially not a single massive cluster but a cluster of numerous small “sub-clusters” that are formed from gas with smaller $`M_\mathrm{J}`$ (significantly smaller than $`10^510^6`$ $`M_{}`$). These clusters are born in the very centers of colliding clouds, and may finally form a single massive cluster. Previous theoretical works argue that the very low $`M_\mathrm{J}`$ (thus very small cluster mass) in the shocked gas layer is a serious problem for the cloud-cloud collision model of star cluster formation (e.g., Kumai et al. 1993). The present numerical results suggest that this Jeans mass problem may not be so important. The incipient sub-clusters may quickly merge with one another to form a single massive cluster due to their compact distribution in the central regions of the colliding clouds.

We may speculate that the long-term evolution of binary clusters, of which a detailed discussion is beyond the scope of this paper, may depend on whether or not the remaining gas is quickly removed during/after the cloud collisions. Such gas removal likely occurs due to the thermal and dynamical effects of young OB stars and type II supernovae. We confirm that if the remaining gas is not removed from the remnant of the cloud-cloud collision in our standard model, the binary cluster eventually merges to form a single cluster because of efficient dynamical friction between the cluster and the low-density gas. In the models with $`b/2r_{\mathrm{cl}}`$ = 0.25, the developed binary clusters coalesces into a single cluster within 0.2 Gyr, for $`V_\mathrm{r}`$ $`<`$ 27 km s-1.

The parameter dependencies are summarized as follows (see Figure 4). Firstly, there is an optimal range of $`V_\mathrm{r}`$ ($`1050`$ km s-1) for star cluster formation. The star cluster formation efficiency rapidly drops as $`V_\mathrm{r}`$ becomes smaller than a threshold value ($``$ 6 km s-1). This occurs because of the much weaker compression and less efficient shock dissipation of the colliding gas. The star formation efficiency becomes very small for models with large $`V_\mathrm{r}`$ ($`>`$ $`50`$ km s-1), since the merging of two clouds does occur at all in these models. Secondly, the star formation efficiency is likely to be higher for models with a smaller impact parameter $`b/2r_{\mathrm{cl}}`$ (in particular, for $`V_\mathrm{r}`$ $`<`$ $`20`$ km s-1). Thirdly, a single cluster rather than a binary cluster is likely to be formed just after the cloud collision in models with smaller $`b/2r_{\mathrm{cl}}`$. Finally, regardless of model parameters, the star clusters possess flattened shapes just after their formation. Our derived higher star formation efficiency can be responsible for the bound cluster formation after gas removal (e.g., Geyer & Burkert 2001).

The present study suggests that the formation of star clusters in colliding gas clouds with $`V_\mathrm{r}`$ $`>`$ 60 km s-1 is much less likely for the adopted parameter range of $`b/2r_{\mathrm{cl}}`$. As is shown in Figure 5, the two colliding clouds with $`b/2r_{\mathrm{cl}}`$ = 0.25 and $`V_\mathrm{r}`$ = 67 km s-1 do not show any star-formation in their central regions for 43 Myr evolution. This is essentially a result of the oblique collision of two clouds with larger $`V_\mathrm{r}`$ ($`>`$ 60 km s-1), they simply graze each other and soon become well separated without forming strongly shocked and compressed gaseous regions conducive to star formation. Therefore it appears that the formation of a thick gaseous layer (with $`M_\mathrm{J}`$ of $`10^510^6`$ $`M_{}`$) in colliding clouds with very large $`V_\mathrm{r}`$ as proposed by Kumai et al. 1993 is unlikely. Our results imply that the formation of star clusters in colliding clouds is more complex than suggested by the previous analytical works of Kumai et al. (1993) and Fujimoto & Kumai (1997).

## 4 Discussions

### 4.1 Comparison with previous works

The dependencies of star cluster formation efficiency on $`V_\mathrm{r}`$ (or $`V_\mathrm{r}/c_\mathrm{s}`$) strongly suggest that the original proposal by Kumai et al. (1993) on collisional cluster formation should be significantly modified. Kumai et al. (1993) adopted the following two assumptions to investigate the formation processes of star clusters in the LMC/SMC system: (1) The relative velocity of two interacting clouds is likely to be more than 50-100 km s-1 and (2) Colliding two gas clouds with such a large relative velocity can form a compressed (shocked) gas layer, where star formation can proceed. Based on these two assumptions, they claimed that the observed differences in cluster formation efficiency between the Galaxy and the LMC/SMC system is due to large-scale random motions with velocities in the Magellanic Clouds, which is not seen in the present-day Galaxy.

We find that the most likely relative velocities ($`V_\mathrm{r}`$) of colliding clouds in the interacting LMC/SMC system is not 50-100 km s-1 but rather $`2540`$ km s-1 for a reasonable set of parameters for the gas clouds. We have also found that only a minor fraction of colliding gas clouds have $`V_\mathrm{r}`$ more than 100 km s-1. The principal reason for these lower velocities than those expected by Kumai et al. (1993) is that the mass ratio of the SMC to the LMC is small enough ($``$ 0.1) that the tidal interaction between the two can not strongly disturb the LMC gas clouds. These results imply that (1) the above first assumption adopted by Kumai et al. (1993) is invalid and (2) if star clusters in LMC/SMC are formed from colliding clouds with $`V_\mathrm{r}`$ more than 50-100 km s-1, these clusters will be a minority among the young star clusters.

Our simulations have showed that if $`V_\mathrm{r}`$ $`>`$ 60 km s-1, the star formation efficiency in colliding clouds becomes very small (therefore any bound star clusters are much less likely to be formed). This is because models with $`V_\mathrm{r}`$ $`>`$ 60 km s-1, result in colliding clouds which just graze with each other and then become well separated without forming a strongly shocked gas layer. Unless the collision of such gas clouds is close to a head-on, the clouds cannot form a bound star cluster. The optimal range of $`V_\mathrm{r}`$ in our simulations suggests that the formation of young star clusters in the LMC/SMC system (but not in the Galaxy) are a result the enhanced rates of cloud-cloud collisions with more moderate relative velocities ($`V_\mathrm{r}`$ $``$ $`1050`$ km s-1). We here stress that this optimal value is true for the clouds with the adopted size-mass relation Larson (1981), masses and sizes in the present study.

Kumai et al. (1993) and Fujimoto & Kumai (1997) found that if gas clouds make a head-on collision with $`V_\mathrm{r}`$ of $``$ 100 km s-1, the Jeans mass ($`M_\mathrm{J}`$) of the compressed thin gas layer becomes about $`10^510^6`$ $`M_{}`$, depending on the sound velocity of the gas layer. The origin of their derived $`M_\mathrm{J}`$ for star cluster formation is their larger adopted value of $`V_\mathrm{r}`$. The present study has demonstrated that star clusters are less likely to be formed in colliding clouds with $`V_\mathrm{r}`$ $``$ 100 km s-1 (also $`M_\mathrm{J}<<10^510^6`$ $`M_{}`$). Therefore, a single star cluster with mass of $`10^510^6`$ $`M_{}`$ is less likely to form from gravitationally unstable gas with $`M_\mathrm{J}`$ of $`10^510^6`$ $`M_{}`$ in colliding clouds with $`V_\mathrm{r}`$ $``$ 100 km s-1.

Our simulations show that star formation starts from the very center of colliding clouds, where $`M_\mathrm{J}`$ is estimated to be $`10^3`$ $`M_{}`$ for our adopted isothermal equation of state. Stars can continue to form in the high-density gas at the cloud center (with $`M_\mathrm{J}`$ of $`10^3`$ $`M_{}`$) such that the central region of each colliding cloud may be regarded as “a cluster of sub-clusters” with the masses of $`10^3`$ $`M_{}`$. These “sub-clusters” are all located in the very center of the colliding cloud, and consequently may quickly merge with one another to form a single massive cluster. Thus, massive clusters ($`10^510^6`$ $`M_{}`$) may result from the merging of numerous small sub-clusters in the cloud centers since $`M_\mathrm{J}`$ of the gas is less than $`10^510^6`$ $`M_{}`$.

Since the resolution of the present simulation is at most $``$ $`10^2`$ $`M_{}`$ in mass and $``$ 1 pc in size for models with $`M_{\mathrm{cl}}`$ = $`10^6`$ $`M_{}`$, we can not investigate the details of the merging process of these sub-clusters. It is necessary to represent such a small cluster not as a single stellar particle, but as a collection of stellar particles (with the total number of $`10^210^3`$), to allow a rigorous investigation of the structural and kinematical properties of the remnants of multiple mergers between numerous small clusters. We leave to a future numerical study, with the total particle numbers of more than $`10^6`$, to confirm (1) whether such small clusters are first formed in the central regions of colliding clouds and (2) how these clusters can merge with one another to form a single massive cluster.

### 4.2 Origin of the observed mass-size relation of star clusters

We have found that the star formation efficiency in colliding gas clouds becomes smaller in the model with larger $`V_\mathrm{r}`$ (or $`V_\mathrm{r}/c_\mathrm{s}`$) for $`V_\mathrm{r}`$ $`>`$ $`30`$ km s-1. This result provides a new clue to the origin of the observed mass-size relation of star and globular clusters. There is observational evidence suggesting only a weak correlation between the mass ($`M_{\mathrm{st}}`$) and size ($`R_{\mathrm{st}}`$) of young star clusters and globular clusters (e.g., $`R_{\mathrm{st}}M_{\mathrm{st}}^{}{}_{}{}^{0.1\pm 0.1}`$ for young clusters; Zepf et al. 1999). Ashman & Zepf (2001) pointed out that if star clusters and globular clusters are formed from molecular gas clouds with the observed size-mass relation of $`r_{\mathrm{cl}}M_{\mathrm{cl}}^{}{}_{}{}^{0.5}`$ (Larson 1981), then the star formation efficiency should be lower in smaller gas clouds to reproduce the observed mass-radius relation. The sound velocity and gas temperature are smaller in smaller self-gravitating gas clouds (Larson 1981). However, $`V_\mathrm{r}`$ is controlled by global galactic dynamics and therefore is independent of gas cloud mass. Our results suggest that star formation efficiency is lower for smaller gas clouds due to the larger $`V_\mathrm{r}/c_\mathrm{s}`$ (for $`V_\mathrm{r}`$ $`>`$ 30 km s-1). We note that the origin of the observed scaling relation of star clusters could be closely associated with the star formation efficiency dependent on $`V_\mathrm{r}/c_\mathrm{s}`$ in colliding clouds.

Following the simple analytic argument by Ashman & Zepf (2001), we can discuss this point in a more quantitative manner. First we define $`ϵ`$ to be the star formation efficiency in colliding gas clouds;

$$ϵ=\frac{M_{\mathrm{st}}}{M_{\mathrm{cl}}}.$$
(1)
We then assume that the size of a cluster depends upon the star formation efficiency of the cluster’s progenitor cloud such that (e.g., Hills 1980);

$$\frac{R_{\mathrm{st}}}{r_{\mathrm{cl}}}ϵ^1.$$
(2)
Based on our simulations, we can write the dependence of $`ϵ`$ on $`V_\mathrm{r}/c_\mathrm{s}`$ as follows;

$$ϵ(\frac{V_\mathrm{r}}{c_\mathrm{s}})^\alpha .$$
(3)
Therefore, the dependence of $`R_{\mathrm{st}}`$ on $`M_{\mathrm{cl}}`$ is:

$$R_{\mathrm{st}}r_{\mathrm{cl}}ϵ^1M_{\mathrm{cl}}^{}{}_{}{}^{1/2+\alpha /4}M_{\mathrm{st}}^{}{}_{}{}^{\frac{2+\alpha }{4\alpha }}$$
(4)
Here we assume that (1) $`r_{\mathrm{cl}}M_{\mathrm{cl}}^{}{}_{}{}^{0.5}`$ (Larson 1981), (2) $`c_\mathrm{s}M_{\mathrm{cl}}^{}{}_{}{}^{0.25}`$ predicted from the virial theorem and Larson’s relation above, and (3) $`V_\mathrm{r}`$ does not depend on initial cloud mass. It is clear from the above equation that $`\alpha `$ should be $``$ $`2`$ to explain the apparent lack of a mass-radius relation in young clusters (Zepf et al. 1999). Although our simulations suggest that $`\alpha `$ takes a negative value, they do not enable us to derive a robust value (or range) for $`\alpha `$ because of our relatively small parameter space. We plan to estimate $`\alpha `$ in more sophisticated future simulations with a much wider parameter space of colliding gas clouds.

### 4.3 Formation of highly flattened star clusters in LMC/SMC

Origin of the flattening of LMC star clusters have been discussed in several authors (e.g., Frenk & Fall 1982; Kontizas et al. 1989). de Oliveira et al. (2000) investigated star cluster encounters in their purely collisionless simulations, and found that binary clusters can merge with each other to form a single cluster with higher ellipticity. They also showed that if the mass ratio of two merging clusters with orbital eccentricities of $`0.60.9`$ (i.e., a bound orbit) is close to 0.1 (“minor merging”), the merger remnant shows an ellipticity consistent with observations of LMC clusters. The present study has demonstrated that binary clusters can be formed from cloud-cloud collisions. Our study, and that of de Oliveira et al. (2000), suggest that the observed globular and populous clusters with higher ellipticity in the Magellanic Clouds may originate from the merging of binary clusters formed from cloud-cloud collisions. These two studies also imply that the difference in the shapes of clusters between the Galactic halo/disk globular clusters and the young Magellanic Cloud clusters may be due to the fact that only LMC/SMC clusters have experienced the past merging of star clusters.

It is, however, unclear whether binary clusters with the mass ratio of $``$ 0.1 can actually be formed in colliding gas clouds. All of our models involve “major mergers” of two clouds such that the incipient star clusters in the centers of the two clouds have similar masses. Also, we have only investigated the collision of two gas clouds with identical radial density profiles. This does not allow us to predict the final mass ratio of two clusters formed from unequal-density mergers of gas clouds. We need to investigate extensively a set of numerical simulations with a wider range of parameters such as the mass ratio and the density profiles of colliding two clouds in order to confirm whether a binary cluster formed from a cloud-cloud collision has the mass ratio of $`sim`$ 0.1. Our future more sophisticated simulations including chemical evolution, magnetic fields, dynamical evolution of hierarchical/fractal structures within a cloud, and feedback effects from massive stars and supernovae will address the origin of flattened shapes of LMC/SMC clusters in a more quantitative way.

## 5 Conclusions

We have used two different sets of numerical simulations to understand the origin of the physical properties of star clusters in the Magellanic Clouds. Although the present model of star cluster formation in colliding two clouds is simplified in a number of areas, we have revealed some essential aspects of star cluster formation in colliding gas clouds.

The main conclusions are summarized as follows.

(1) An oblique collision between two identical clouds can be significantly enhanced during the tidal interaction between the LMC and the SMC. Cloud-cloud collisions with radial velocities ($`V_\mathrm{r}`$) of $`2540`$ km s-1 are most common in our models. These results imply that the origin of young star clusters in the LMC/SMC systems may stem from the enhanced collision rates of gas clouds with moderately large relative velocities. Our results also suggest that the original proposal by Kumai et al. (1993) on collisional cluster formation should be modified.

(2) Compact, bound star clusters can be formed in the centers of colliding gas clouds as a result of strong gas shocks, compression, and dissipation during the collision. The initial impact parameter of two colliding clouds determines whether the star clusters form a single cluster, a binary cluster, or two isolated star clusters. For example, for a smaller impact parameter, the incipient clusters soon merge with each other to form a single, more massive cluster.

(3) Star formation efficiency in colliding clouds depends on the initial ratio of the relative velocity of the clouds to the sound speed of the gas ($`V_\mathrm{r}/c_\mathrm{s}`$). This dependency is in the sense that the star formation efficiency is lower for models with larger $`V_\mathrm{r}/c_\mathrm{s}`$. The derived dependence on $`V_\mathrm{r}/c_\mathrm{s}`$ provides a new clue to the origin of the observed mass-size relation of young star clusters.

We are grateful to the anonymous referee for valuable comments, which contribute to improve the present paper. K.B. and W.J.C. acknowledge the Large Australian Research Council (ARC). MB acknowledges assistance from the Swinburne Research and Development Grant System.