text:
# Star clusters and galactic chemodynamics: Implosive formation of super star clusters.

## 1 Introduction

Star clusters (hereafter SC) such as open, globular, and super-star clusters are observed to be relatively “minor” components (in mass) in luminous galaxies (e.g., the mass ratio of the Galactic GCs to the Galactic luminous disk mass is an order of 0.1 %). Therefore, it is generally considered that dynamical evolution of star clusters (e.g., merging or tidal destruction of star clusters) does not influence significantly the evolution of luminous galaxies. However, recent numerical simulations on dynamical evolution of SCs have strongly suggested that this SC dynamics is very important for the formation of galactic nuclei in dwarf irregulars/ellipticals and that of ultra-compact dwarf galaxies (Drinkwater et al. 2003), because (1) time scale of dynamical friction of SCs is relatively short ($``$ $`10^9`$ yr) and (2) cross section of SC merging is also large owing to the larger $`\mathrm{S}_\mathrm{N}`$ of dwarf galaxies (Bekki et al. 2003). Furthermore recent HST observations have revealed a significant number of young super star clusters (hereafter SSC) in dwarf irregular galaxies and suggested that physical properties of these SSCs can be correlated with those of their host galaxies (e.g., Billett et al. 2002). These recent numerical and observational results lead the author to investigate (1) how SSCs can be formed in galaxies, (2) how these SSCs dynamically evolve, and (3) whether or not there can be strong physical relationship between structural and chemical properties of SSCs and those of their host galaxies. In this paper, we discuss these problems based on numerical simulations of SSC formation in galaxies. We here focus on SSC formation in major mergers, because recent HST observations discovered many SSCs in interacting/merging galaxies (e.g., Whitmore et al. 1993).

## 2 A new two-fold numerical model

Since it is very hard to resolve SSC formation sites (less than pc) and investigate the formation processes in detail in large-scale (10-100kpc) galaxy simulations, we adopt the following novel two-fold method for the numerical investigation of SSC formation. Firstly, we investigate time evolution of pressure, temperature, and density of warm interstellar medium (ISM) with the initial temperature of $`10^4`$ K and the density of 1 atom $`\mathrm{cm}^3`$, based on galaxy-scale (larger-scale) SPH simulations on galaxy mergers with total particle number ($`N`$) of $`10^5`$. Secondly, we investigate dynamical evolution and star formation histories of GMCs embedded by hot (warm), high pressure ISM, based on GMC-scale ($`<`$ 100pc) SPH simulations ($`N4\times 10^4`$). Primary focus here is to investigate how the high pressure ISM (expected for galaxy mergers) trigger collapse of GMCs and consequently form SSCs within the GMCs. In the second simulations, the data of the ISM surrounding GMCs come from the first simulations: The present simulations are fully self-consistent and mimic (or corresponds to) the simulations with $`N10^9`$.

### 2.1 Galaxy-scale

Our numerical methods for modeling chemodynamical evolution of galaxy mergers have already been described by Bekki & Couch (2001) and the details of the adopted TREESPH are given by Bekki (1995). We accordingly give only a brief review of them here. We construct models of galaxy mergers between gas-rich disks by using the model of Fall-Efstathiou (1980). The total mass and the size of a progenitor exponential disk are $`M_\mathrm{d}`$ and $`R_\mathrm{d}`$, respectively. From now on, all the masses and lengths are measured in units of $`M_\mathrm{d}`$ and $`R_\mathrm{d}`$, respectively, unless otherwise specified. Velocity and time are measured in units of $`v`$ = $`\left(GM_\mathrm{d}/R_\mathrm{d}\right)^{1/2}`$ and $`t_{\mathrm{dyn}}`$ = $`\left(R_\mathrm{d}^3/GM_\mathrm{d}\right)^{1/2}`$, respectively, where $`G`$ is the gravitational constant and assumed to be 1.0 in the present study. If we adopt $`M_\mathrm{d}`$ = 6.0 $`\times `$ $`10^{10}`$ $`\mathrm{M}_{}`$ and $`R_\mathrm{d}`$ = 17.5 kpc as fiducial values, then $`v`$ = 1.21 $`\times `$ $`10^2`$ km $`\mathrm{s}^1`$ and $`t_{\mathrm{dyn}}`$ = 1.41 $`\times `$ $`10^8`$ yr, respectively. The dark-to-disk halo mass ratio and the star-to-gas mass ratio are set equal to 4.0 and 9.0, respectively. Bulge component is not included in the present study. An isothermal equation of state is used for the gas with a temperature of $`7.3\times 10^3`$ K (corresponding to a sound speed of 10 km $`\mathrm{s}^1`$). We present the results mostly for merger models with nearly prograde-prograde orbital configurations, parabolic encounters, and a pericentric distance of 0.5 $`R_\mathrm{d}`$.

### 2.2 GMC-scale

We numerically investigate the hydrodynamical effects of the high pressure ISM on a self-gravitating molecular gas cloud in major mergers. The cloud is represented by 20,000 SPH particles and the initial cloud mass ($`M_{\mathrm{cl}}`$) and size ($`r_{\mathrm{cl}}`$) are set to be $`10^6M_{}`$ and 97 pc, respectively, which are consistent with the mass-size relation observed by Larson (1981). The cloud is assumed to have an isothermal radial density profile with $`\rho \left(r\right)1/\left(r+a\right)^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 be 4 km s-1 (consistent with the prediction from the virial theorem) for models with $`M_{\mathrm{cl}}=10^6M_{}`$. Each SPH particle is also subject to the hydrodynamical force of the ISM whose strength depends on the parameters of the ISM. The ISM is represented by $``$ 20000 SPH particles and and an isothermal equation of state with pressure $`p`$ (or sound velocity of $`c_\mathrm{h}`$), and density, $`\rho `$. The ISM particles are uniformly distributed within a box with size of $`6R_{\mathrm{cl}}`$. A gas particle in a cloud is converted into a collisionless stellar particle if (1) the local dynamical time scale \[corresponding to $`\left(4\pi G\rho _i\right)^{0.5}`$, where $`G`$ and $`\rho _i`$ are the gravitational constant and the density of the gas particle, respectively\] is shorter than the sound crossing time (corresponding to $`h_i/c_\mathrm{s}`$, where $`h_i`$ is the smoothing length of the gas), and (2) the gaseous flow is converging. This method thus mimics star formation due to Jeans instability in gas clouds. In the model without ISM effects, star formation does not occur at all in this star formation method. We mainly show the result of the “fiducial model” with$`\rho `$ = 1 atom cm-3 and $`P`$ = $`10^5`$ $`\mathrm{k}_\mathrm{B}`$ K cm-3, because this model shows the typical behavior of ISM-induced star formation in a gas cloud.

## 3 Result I. SSC formation sites.

Figure 1 clearly demonstrates how the gaseous pressure of the interstellar medium becomes dramatically higher ($`10^4`$-$`10^5`$ $`k_\mathrm{B}`$) in some regions of the major merger model at $`T`$ = 0.28 Gyr (starburst phase with the star formation rate of $``$ 75 $`M_{}`$ $`\mathrm{yr}^1`$) in comparison with the initial disks. Furthermore, gaseous pressure for some fraction of the gas particles (in particular, in the central regions) exceeds the threshold value of $`2.0\times 10^5k_\mathrm{B}`$ for SSC formation, essentially because rapid transfer of gas to the central region of the merger and efficient gaseous shock dissipation increase the gaseous density and pressure greatly. These results clearly explain the reason why SSCs can be formed in starbursting galaxies and imply that the formation efficiency of SSCs are much higher in merging galaxies with massive starbursts than in “normal” disk galaxies. Figure 1 also shows that SSCs are formed not only in the central region, where gas fueling is more efficient, but also in the outer tidal tails and “bridges” among the two cores, where the tidal force of the merger forms thin and high-density gaseous layers owing to gaseous shock dissipation (More details are given in Bekki & Couch 2001).

## 4 Result II. Implosive SSC formation.

Figures 2 and 3 describe how a burst of star formation within a gas cloud can be triggered by the external gas pressure of the ISM in the fiducial model. In this model with moderately strong ISM pressure, the high pressure of the ISM can continue to strongly compress the cloud without losing a significant amount of gas from the cloud. As the strong compression proceeds, the internal density/pressure of the cloud can rise significantly. However, the self gravitational force of the cloud also becomes stronger because the cloud becomes progressively more compact during the ISM’s compression. Therefore, the internal gaseous pressure of the cloud alone becomes unable to support itself against the combined effect of the external pressure from the ISM and the stronger self-gravitational force. As a result of this, the cloud’s collapse initially induced by the high external pressure from the ISM can continue in a runaway manner.

Due to the rapid, dissipative collapse, the gaseous density of the cloud dramatically rises and consequently star formation begins in the central regions of the cloud. The star formation rate increase significantly to 1.5 $`M_{}`$ yr-1 (8 Myr after the start of the cloud’s collapse). About 80 % of the gas is converted into stars within 14 Myr to form a stellar system. Because of the “implosive” formation of stars from strongly compressed gas, the developed stellar system is strongly self-gravitating and compact. This result implies that high external pressure from the ISM is likely to trigger the formation of bound, compact star clusters rather than unbound, diffuse field stars. The final mass distribution of this model is similar to the King profile, but too compact (with the effective radius of $``$ 2 pc) to be consistent with the observed effective radius of typical globular clusters (See Figure. 4).

## 5 Conclusions.

We first demonstrated that SSCs can be formed from GMCs “in an implosive manner” for major mergers: Strong external pressure of ISM in mergers can trigger the collapse of GMCs and consequently cause the rapid star formation within the GMCs. The star formation efficiency in the GMCs under strong external pressure is as high as $`6080`$ % (in total cloud mass) so that the developed SSCs can be considered to be bound star clusters (e.g., Hills 1980). The present study suggests that unbound star clusters that are doomed to be destroyed by some external processes (e.g., galactic tidal field and interaction with GMC) can be formed from GMCs embedded by low (normal) pressure ISM whereas SSCs can be formed from GMCs embedded by very high pressure ISM. The stars initially within unbound star clusters could soon be identified as “field stars” after the destruction of the clusters. The time scale of SSC formation is very short (a few times 106 yr, which is shorter than the lifetime of very massive stars, which are responsible for the type II supernovae) so that the energy feedback from supernovae can be much less likely to influence the SSC formation. Although the “implosive formation scenario” described here can provide a clue to the origin of GCs in general, we have not yet reproduced the observed Fundamental Plane of GCs (in M31 and the Galaxy) and mass-size relation of GCs. Therefore, it is our future studies to check whether the implosive scenario can explain self-consistently the observed scaling relations of GCs and chemical properties of GCs.

## Acknowledgments

## References

Bekki, K. 1997, ApJ, 483, 608

Bekki, K., Couch, W. J. 2001, ApJL, 557, 19

Bekki, K., Couch, W. J., Drinkwater, M. J., Shioya, Y. 2003, MNRAS, in press.

Billett, O., Hunter, D. A., Elmegreen, B. G., 2002, AJ, 123, 1454

Drinkwater, M. J., Gregg, M. D., Hilker, M., Bekki, K., Couch, W. J., Ferguson, H. C., Jones, J. B., Phillipps, S. 2003, Nature, 423, 519

Fall, S. M., & Efstathiou, G. 1980, MNRAS, 193, 189

Hills, J. G. 1980, ApJ, 235, 986

Larson, R. B. 1981, MNRAS, 194, 809

Whitmore, B. C., Schweizer, F. Leitherer, C. Borne, K., Robert, C. 1993, AJ, 106, 1354


format: text

text:
# Drifting diffusion on a circle as continuous limit of a multiurn Ehrenfest model

## Abstract

We study the continuous limit of a multibox Erhenfest urn model proposed before by the authors. The evolution of the resulting continuous system is governed by a differential equation, which describes a diffusion process on a circle with a nonzero drifting velocity. The short time behavior of this diffusion process is obtained directly by solving the equation, while the long time behavior is derived using the Poisson summation formula. They reproduce the previous results in the large $`M`$ (number of boxes) limit. We also discuss the connection between this diffusion equation and the Schr$`\ddot{\mathrm{o}}`$dinger equation of some quantum mechanical problems.

In a previous study KL  we proposed a generalized Ehrenfest urn model Erf  with $`N`$ balls and $`M`$ urns that arranged periodically along a circle. The evolution of the system is governed by a directed stochastic operation. Using the standard matrix diagonalization procedures together with a multi-variable generating function method, we have solved the problem completely. We found that for a generic $`M>2`$ case the average number of balls in a certain urn oscillates several times before it reaches a stationary value. We also obtained the Poincaré cycle Huang , i.e., the average time interval required for the system to return to its initial configuration. The result is simply given by $`M^N`$, which indicates that the fundamental assumption of statistical mechanics holds in this system. Taking $`M=2`$, our model reproduces all the results of the original Erhenfest urn model Erf .

In this paper, we further study the continuous limit (the large $`M`$ and $`N`$ limit) of the proposed multiurn model. We show that by defining a density function $`\rho `$ as the continuous limit of the fraction $`f_i=m_i/N`$, i.e., the average number of balls in the $`i`$th urn devided by $`N`$, the continuous limit of the model exists if we also define the drifting velocity and diffusion constant appropriately. The evolution of $`\rho `$ in spacetime is then governed by a differential equation, which can be solved under proper initial condition and boundary conditions. The results obtained in this paper are in agreement with those obtained before by the standard matrix diagonalization method. Since for even a generic $`M`$-urn and $`N`$-ball case the Poincar’e cycle $`M^N`$ is too huge to be experienced, the evolution of the system can in practice be treated as unrepeatable, thus the average quantities considered here become more important than those of microstate details.

We start from the Eq. (4) of Ref.KL :

$$m_i_s=\left(1\frac{1}{N}\right)m_i_{s1}+\frac{1}{N}m_{i1}_{s1},$$
(1)
which can be rewritten as

$$f_i(s)f_i(s1)=\frac{1}{N}\left[f_i(s1)f_{i1}(s1)\right],$$
(2)
where $`f_i(s)m_i_s/N`$, $`N`$ is the total number of the balls, and $`m_i_s`$ denotes the number of balls in the $`i`$th urn after $`s`$ steps. Adding $`[f_{i+1}(s1)f_{i1}(s1)]/2N`$ to both sides of Eq. (2), we get

$`{\displaystyle \frac{f_i(s)f_i(s1)}{\mathrm{\Delta }t}}+{\displaystyle \frac{\mathrm{\Delta }x}{N\mathrm{\Delta }t}}\left[{\displaystyle \frac{f_{i+1}(s1)f_{i1}(s1)}{2\mathrm{\Delta }x}}\right]`$
$`={\displaystyle \frac{(\mathrm{\Delta }x)^2}{2N\mathrm{\Delta }t}}\left[{\displaystyle \frac{f_{i+1}(s1)2f_i(s1)+f_{i1}(s1)}{(\mathrm{\Delta }x)^2}}\right],`$ (3)
where $`\mathrm{\Delta }t`$ represents the time interval in one step, and $`\mathrm{\Delta }x`$ stands for the center-center distance between two neighboring urns. Taking the continuous limit, we obtain

$$\frac{\rho }{t}+v\frac{\rho }{x}=D\frac{^2\rho }{x^2},$$
(4)
where we have used the substitutions:

$$f_i(s)\rho (x,t),\frac{\mathrm{\Delta }x}{N\mathrm{\Delta }t}v,\frac{(\mathrm{\Delta }x)^2}{2N\mathrm{\Delta }t}D.$$
(5)
It is clear that Eq. (4) is a diffusion equation. Since the model is defined on a circle, we replace $`x`$ by $`\varphi `$, $`v`$ by $`\omega `$, $`\mathrm{\Delta }x`$ by $`\theta `$, and the diffusion equation becomes

$$\frac{\rho }{t}+\omega \frac{\rho }{\varphi }=D\frac{^2\rho }{\varphi ^2}.$$
(6)
Before further exploring Eq.(6), here we give a simple and general derivation of the diffusion equation. Note that the conservation of probability implies

$$\frac{\rho }{t}=𝐉,$$
(7)
where $`\rho (𝐫,t)`$ is the probability density and $`𝐉(𝐫,t)`$ is the probability current density. Now, the probability current can be written as the sum of two terms, one for the “diffusion part”, and the other for the “drifting part” of the probability carriers (the balls). That is

$$𝐉=D\rho +\rho 𝐯,$$
(8)
where $`D`$ is the diffusion constant and $`𝐯`$ is the drifting velocity caused by some pumping force.

Substitute (8) into (7), we obtain

$$\frac{\rho }{t}=D^2\rho (\rho 𝐯).$$
(9)
We further assume that $`𝐯=0`$ (incompressible fluid; one special case is that $`𝐯=`$ constant), then we have

$$\frac{\rho }{t}=D^2\rho 𝐯\rho ,$$
(10)
which is the desired diffusion equation and has the same form as Eq. (4) and (6).

On a straight line, the above equation becomes Eq. (4), and we adopt the boundary condition

$$\rho (\mathrm{},t)=\rho (\mathrm{},t)=0.$$
(11)
On a circle, Eq. (10) becomes Eq. (6), with boundary condition

$$\rho (\varphi ,t)=\rho (\varphi +2\pi ,t).$$
(12)
Now we find the solutions $`\rho `$ for the 1D diffusion equations on a straight line (4) and on a circle (6), respectively. Assuming the initially condition:

$$\rho (x,0)=\delta (x),$$
(13)
the solution on a line can be obtained by Fourier transform method MW :

$$\rho (x,t)=\frac{1}{\sqrt{4\pi Dt}}\mathrm{exp}\left[\frac{(xvt)^2}{4Dt}\right].$$
(14)
Similarly, for the circle problem, given the initial condition

$$\rho (\varphi ,0)=\delta (\varphi ),$$
(15)
we obtain

$$\rho (\varphi ,t)=\frac{1}{\sqrt{4\pi Dt}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left[\frac{(\varphi \omega t+2n\pi )^2}{4Dt}\right].$$
(16)
In deriving Eq. (16), we have used the identity

$$_{\mathrm{}}^{\mathrm{}}f(x)𝑑x=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}_0^{2\pi }f(x+2n\pi )𝑑x$$
(17)
for a localized function $`f(x)`$, and we have treated the “circle problem” as an “infinite-folded line problem”.

Furthermore, the “center of mass” now is written

$`\mathrm{COM}`$ $`=`$ $`{\displaystyle _0^{2\pi }}𝑑\varphi \rho (\varphi ,t)\mathrm{exp}(i\varphi )`$ (18)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{4\pi Dt}}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _0^{2\pi }}𝑑\varphi e^{\frac{(\varphi +2n\pi \omega t)^2}{4Dt}+i(\varphi +2n\pi )}`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{4\pi Dt}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\varphi e^{\frac{(\varphi \omega t)^2}{4Dt}+i\varphi }`$
$`=`$ $`{\displaystyle \frac{e^{i\omega t}}{\sqrt{4\pi Dt}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\varphi e^{\frac{\varphi ^2}{4Dt}+i\varphi }`$
$`=`$ $`\mathrm{exp}(Dt+i\omega t),`$
which is equivalent to Eq.(32) of Ref. KL  if we define

$$Dt=\frac{\theta ^2}{2}\tau =\frac{2\pi ^2}{M^2}\tau ,\omega t=\theta \tau =\frac{2\pi }{M}\tau .$$
(19)
Here $`\tau `$ and $`\theta `$ are defined as

$$\tau \frac{t}{N\mathrm{\Delta }t}=\frac{s}{N},\theta =\mathrm{\Delta }x=\frac{2\pi }{M}.$$
(20)
Now we compare the results with those in Ref. KL . Fig. (1) shows the results from Eq. (16), (19), and (20) for the cases $`M=30`$ and $`M=60`$. As one can see, they indeed reproduce the results of Ref. KL  in the large $`M`$ limit. The parameter $`N`$ does not appear here because the motion of each particle is independent in our model.

Note that although for a small $`t`$ the expression (16) is good enough to be a fast convergent series, however, when $`t`$ becomes large, Eq. (16) converges slowly. In this situation we use a more accurate expression for $`\rho `$:

$`\rho (\varphi ,t)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{n^2Dt+in(\varphi \omega t)}`$ (21)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}+{\displaystyle \frac{1}{\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}e^{n^2Dt}\mathrm{cos}[n(\varphi \omega t)],`$
which can be derived from (16) using the Poisson summation formula: CH

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}f(na)=\frac{2\pi }{a}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}g\left(\frac{2n\pi }{a}\right).$$
(22)
Here $`f(x)`$ is a localized function, and

$$g(k)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}f(x)e^{ikx}𝑑x$$
(23)
is its Fourier transform.

We now consider some solvable generalizations of Eq. (4) and (6). Note that the ratio between $`D`$ and $`\omega `$ in our model is fixed:

$$\frac{D}{\omega }=\frac{\theta }{2}=\frac{\pi }{M}.$$
(24)
To relax this restriction, we modify our urn model by assuming that at each time step the picked ball can have probability $`p`$ to be put into the next urn and probability $`q=1p`$ to be put into the previous urn. Hereafter we call this modified model the $`pq`$-model. The $`pq`$-model is also solvable KL1  by using methods like those used in Ref. KL . The continuous limit of the $`pq`$-model can be derived from the recurrence relation for $`f_i`$:

$`f_i(s)`$ $`=`$ $`\left(1{\displaystyle \frac{1}{N}}\right)f_i(s1)`$ (25)
$`+`$ $`{\displaystyle \frac{p}{N}}f_{i1}(s1)+{\displaystyle \frac{q}{N}}f_{i+1}(s1),`$
where

$$p+q=1,\mathrm{and}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}p,q1.$$
(26)
Adding

$$\frac{(pq)}{2N}[f_{i+1}(s1)f_{i1}(s1)]$$

to both sides of Eq. (25), it becomes

$`{\displaystyle \frac{f_i(s)f_i(s1)}{\mathrm{\Delta }t}}`$
$`+{\displaystyle \frac{2(pq)\mathrm{\Delta }x}{2N\mathrm{\Delta }t}}\left[{\displaystyle \frac{f_{i+1}(s1)f_{i1}(s1)}{2\mathrm{\Delta }x}}\right]`$
$`={\displaystyle \frac{(\mathrm{\Delta }x)^2}{2N\mathrm{\Delta }t}}\left[{\displaystyle \frac{f_{i+1}(s1)2f_i(s1)+f_{i1}(s1)}{(\mathrm{\Delta }x)^2}}\right].`$ (27)
Now, define

$$f_i(s)\rho (x,t),\frac{(pq)\mathrm{\Delta }x}{N\mathrm{\Delta }t}v,\frac{(\mathrm{\Delta }x)^2}{2N\mathrm{\Delta }t}D,$$
(28)
then we get a continuous equation of the form (4), without the restriction (24). One special case is $`p=q=1/2`$, which has a zero drifting velocity, and the evolution of the system is governed by pure diffusion process — the random walk.

For another generalization we assume that the drifting velocity $`v`$ varies with time, that is

$$\frac{\rho }{t}+v(t)\frac{\rho }{x}=D\frac{^2\rho }{x^2}.$$
(29)
Defining $`x(t)`$ as the time integral of $`v(t)`$:

$$x(t)=_0^tv(t^{})𝑑t^{},$$
(30)
and adopting the initial condition (13), then

$$\rho (x,t)=\frac{1}{\sqrt{4\pi Dt}}\mathrm{exp}\left[\frac{\left(xx(t)\right)^2}{4Dt}\right].$$
(31)
Similarly, for the diffusion equation on a circle with a time-dependent $`\omega (t)`$ and initial condition (15):

$$\frac{\rho }{t}+\omega (t)\frac{\rho }{\varphi }=D\frac{^2\rho }{\varphi ^2},$$
(32)
and the solution is

$$\rho (\varphi ,t)=\frac{1}{\sqrt{4\pi Dt}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left[\frac{(\varphi \varphi (t)+2n\pi )^2}{4Dt}\right].$$
(33)
Here

$$\varphi (t)=_0^t\omega (t^{})𝑑t^{}.$$
(34)
The reason for why $`v`$ and $`\omega `$ can freely vary with time relies on Eq.(28). Recall that in our original multiurn Ehrenfest model or the $`pq`$-model both the time interval between two steps and the distance (angle difference) between two urns are undefined. Thus in deriving the continuous limit of these models we do not have to adopt a constant $`\mathrm{\Delta }t`$ at each step or a fixed $`\mathrm{\Delta }x`$ ($`\mathrm{\Delta }\theta `$) between two neighboring urns. If we relax the restriction in Eq. (28) and modify them to $`\mathrm{\Delta }t_s`$ and $`\mathrm{\Delta }x_i`$ ($`\mathrm{\Delta }\theta _i`$), then the continuous limit of these quantities lead to $`v(t)`$ or $`\omega (t)`$.

It is interesting to note that the solutions for the diffusion equation (10) can be used to find the wave function or Green’s function of some time-dependent quantum mechanical problems Bauer . The main idea is to define a transformation appropriately between the parameters used in the diffusion equation (10) or (4) and (6) and those used in the corresponding Schr$`\ddot{\mathrm{o}}`$dinger equations. For instance, consider a quantum point particle of charge $`q`$ and mass $`m`$ moving under the influence of a vector potential $`𝐀(t)`$ Feynman :

$$i\mathrm{}\frac{\psi }{t}=\frac{1}{2m}\left(i\mathrm{}\frac{q𝐀(t)}{c}\right)^2\psi ,$$
(35)
here we have assumed that $`𝐀(t)`$ is a function of time $`t`$ only. Rewriting Eq. (35) as

$$\frac{\psi }{t}=\frac{i\mathrm{}}{2m}^2\psi +\left(\frac{q𝐀}{mc}\right)\psi \frac{iq^2A^2}{2\mathrm{}mc^2}\psi $$
(36)
and executing the transformation

$$\psi =\mathrm{exp}\left[\frac{1}{2i\mathrm{}m}_0^t\left(\frac{qA(t^{})}{c}\right)^2𝑑t^{}\right]\stackrel{~}{\psi },$$
(37)
we find

$$\frac{\stackrel{~}{\psi }}{t}=\frac{i\mathrm{}}{2m}^2\stackrel{~}{\psi }+\left(\frac{q𝐀}{mc}\right)\stackrel{~}{\psi }.$$
(38)
Comparing Eq. (38) with (10), we find that they can be transformed to each other by the substitution:

$$D\frac{i\mathrm{}}{2m},𝐯\left(\frac{q𝐀}{mc}\right),\rho \stackrel{~}{\psi }.$$
(39)
To be more specific, consider the case that the particle moving on a circle of radius $`1`$. Suppose the circle is lying on the xy-plane and centered at $`(x,y)=(0,0)`$. The vector potential can be chosen as $`𝐀(t)=A(t)\widehat{e}_\varphi `$ and is generated by a time-dependent magnetic flux $`\mathrm{\Phi }(t)`$ tube going through the origin and pointing along the $`z`$-axis

$$𝐀(t)=A(t)\widehat{e}_\varphi =\frac{\mathrm{\Phi }(t)}{2\pi }\widehat{e}_\varphi .$$
(40)
Choosing the initially condition as

$$\stackrel{~}{\psi }(\varphi ,0)=\psi (\varphi ,0)=\delta (\varphi ),$$
(41)
then

$$\psi (\varphi ,t)=\frac{U(t)}{\sqrt{2\pi i\mathrm{}t/m}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left[\frac{(\varphi \varphi (t)+2n\pi )^2}{(2i\mathrm{}t/m)}\right].$$
(42)
Here

$$U(t)=\mathrm{exp}\left[\frac{1}{2i\mathrm{}m}_0^t\left(\frac{qA(t^{})}{c}\right)^2𝑑t^{}\right],$$
(43)
and

$$\varphi (t)=\frac{q}{mc}_0^tA(t^{})𝑑t^{}.$$
(44)
Note that the $`\psi `$ in Eq. (42) is nothing but the Green’s function $`G(\varphi ,\varphi _0;t,t_0)`$ for the quantum particle with $`\varphi _0=t_0=0`$. If $`A(t)=0`$, Eq. (42) gives the well known results HK :

$$G(\varphi ;t)=\frac{1}{\sqrt{2\pi i\mathrm{}t/m}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left[\frac{(\varphi +2n\pi )^2}{(2i\mathrm{}t/m)}\right],$$
(45)
and

$`G(\varphi ;t)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{n^2(i\mathrm{}/2m)t+in\varphi }`$ (46)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}+{\displaystyle \frac{1}{\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}e^{n^2(i\mathrm{}/2m)t}\mathrm{cos}n\varphi `$
for small and large $`t`$, respectively.

In conclusion, we have derived the continuous limit of a mutiurn Ehrenfest model, which is a diffusion equation with a drifting velocity term. Solving the equation gives us the correct time evolution behavior of the ball distribution. A transformation was introduced, which changes the solution of the diffusion equation to the corresponding solution for the problem of a quantum particle moving under the influence of a time-varying magnetic field.

The supports from NSC, NCTS and NCU are thanked. Discussions with Dr. D. H. Lin and Dr. C. S. Tang are also acknowledged.


text:
# Pairing and Condensation of Laughlin Quasiparticles in Fractional Quantum Hall Systems

## I Introduction

In a recent experiment pan , Pan et al. observed the fractional quantum Hall (FQH) effect tsui ; laughlin  at novel filling fractions $`\nu `$ of the lowest Landau level (LL). Among other features, the FQH effect is a macroscopic manifestation of an incompressible character of the many-body ground state formed at a specific filling $`\nu `$. The new spin-polarized FQH states occur at filling factors outside the Jain sequence jain  of composite fermion (CF) states. Some of them, as $`\nu =4/11`$ or 4/13, appear in the Haldane hierarchy haldane  of quasiparticle (QP) condensates, and therefore at first sight it might only be surprising that they have not been observed earlier. However, the “hierarchical” interpretation of these states was questioned hierarchy  because of the specific form of the QP–QP interaction, which leaves their interpretation quite uncertain. Others, such as the $`\nu =3/8`$ or 3/10 states, do not belong to the Haldane hierarchy, and thus the origin of their incompressibility is puzzling in an even more obvious way. In the CF picture, these even-denominator fractions correspond to the half-filled first excited ($`n=1`$) CF LL. While for the electrons half-filling their $`n=1`$ LL (including the double degeneracy of the lowest LL, this amounts to the total electron filling of $`\nu =2+1/2=5/2`$) the FQH effect has been known for quite some time willet  (the incompressible $`\nu =5/2`$ state is known as a Moore-Read state moore ), the similar behavior of electrons and CF’s at this filling is rather unexpected. Pan et al. take their observations as evidence for residual CF–CF interactions, sufficiently strong to cause emergence of new FQH states (that, unlike virtually all observed fractions, cannot be predicted from the noninteracting CF model or standard Haldane hierarchy). However, they ignore the theoretical investigations in which these interactions were studied in detail hierarchy ; sitko ; lee .

In this paper we propose an explanation for these new states involving formation of QP pairs which display Laughlin correlations with one another. First, we explain the connection between the CF model and the QP hierarchy haldane ; hierarchy ; sitko . Second, we recall two simple types of two-body correlations, Laughlin correlations and pairing, that may occur in an interacting system depending on the filling factor $`\nu `$ and on whether $`V()`$ is super- or subharmonic at the relevant range parentage . Then, knowing the QP–QP pseudopotential $`V_{\mathrm{QP}}()`$, we apply the concept of Laughlin condensed states of (bosonic) pairs (used earlier for the electrons in the $`n=1`$ LL to describe such FQH states as $`\nu =5/2`$ or 7/3 fivehalf ) to the particles or holes in a partially filled CF LL, i.e., to Laughlin quasielectrons (QE’s) or quasiholes (QH’s). Finally, we propose the novel hierarchy of FQH states in which the incompressibility results from the condensation of QP pairs (QE2’s or QH2’s) into Laughlin correlated pair states. The series of FQH states derived from the parent $`\nu =1/3`$ state include all novel fractions: $`\nu =5/13`$, 3/8, 4/11, and 6/17 for the QE’s and $`\nu =5/17`$, 3/10, 4/13, and 6/19 for the QH’s.

## II Laughlin QP’s and CF’s

Let us begin with recalling the connection between Laughlin QP’s laughlin  and the CF model (equivalent to the mean-field Chern–Simons transformation) jain . Laughlin QP’s are the actual elementary excitations of a two-dimensional electron liquid filling a fraction of the lowest LL. They have well-defined and known wave functions (and thus also density profiles, size, etc.) and single-particle energy. They carry (fractional) electric charge, counted with respect to the uniformly spread charge of the underlying Laughlin state. The negatively and positively charged QP’s are called quasielectrons (QE’s) and quasiholes (QH’s), respectively. Moving in the external magnetic field, both kinds of QP’s follow cyclotron orbits, conveniently labeled by angular momentum.

Since the uniform-density Laughlin states only occur at a discrete series of filling factors ($`\nu =1/3`$, $`1/5`$, …), and Laughlin QP’s are simply the least-energy excitations able to carry charge in excess of the value corresponding to the nearest Laughlin state, their type and number depends on $`\nu `$. At $`\nu `$ precisely equal to any of the Laughlin values $`\nu _p=(2p+1)^1`$ no QP’s exist in the ground state and they may only appear in form of excitonic, charge-neutral QE–QH pairs. At $`\nu `$ slightly smaller (or larger) than $`\nu _p`$, a number $`N_{\mathrm{QP}}`$ of the “$`p`$–type” QE’s (or QH’s) appear in the Laughlin liquid. The QP number is equal to the difference between the total magnetic flux through the sample and its value corresponding to the nearest “parent” Laughlin state, measured in units of elementary flux, $`N_{\mathrm{QP}}=|\mathrm{\Phi }\mathrm{\Phi }_p|/\varphi _0`$.

For electrons, the filling factor is expressed through the flux $`\mathrm{\Phi }`$ and the electron number $`N`$ as $`\nu =\mathrm{\Phi }/N\varphi _0`$, For the QP’s of the parent $`\nu _p`$ Laughlin state, each carrying charge $`\pm e/(2p+1)`$, the LL degeneracy is reduced by a factor of $`(2p+1)`$ compared to the electron value. Consequently, the QP filling factor is equal to $`\nu _{\mathrm{QP}}=\mathrm{\Phi }/(2p+1)N_{\mathrm{QP}}\varphi _0`$. Note that $`\nu _{\mathrm{QP}}`$ is a linear function of $`\nu `$ near each value of $`\nu _p`$, The situation with $`\nu `$ far away from the nearest $`\nu _p`$, corresponds to a large filling factor for appropriate QP’s whose mutual interactions and filling of higher QP LL’s both become important.

In the mean-field CF picture, the reduction of the QP LL degeneracy is attributed to a reduced magnetic field $`B^{}=B/(2p+1)`$ rather than to a reduced QP charge. In the simplest formulation of the model, the effective field $`B^{}`$ results from the capture of an even number of magnetic flux quanta $`2p\varphi _0`$ by each electron (such bound state is called a CF), and the QE’s and QH’s are pictured as particles and vacancies in the otherwise empty (full) CF LL’s.

The CF picture turns out very useful for the description of many properties of the QP’s (e.g., size of cyclotron orbits or LL degeneracy) or of the FQH systems in general (e.g., the values of $`\nu `$ at which the incompressibility results from a complete filling of the QP/CF LL’s). However, the QP–QP interactions important for the present problem in the CF model arise as a combination of rather complicated two- and three-body gauge interactions between charges and fluxes, and because of this difficulty, are usually neglected. It is therefore important to realize that the QP–QP interaction really is a Coulomb interaction between a pair of charged particles. For example, for two identical QP’s it is repulsive, and at long range it is similar to the electron–electron repulsion, only reduced in magnitude by $`(2p+1)^2`$ because of a smaller QP charge, On the other hand, the exact form of the QP–QP interaction at short range (where it is different from the electron–electron repulsion because of the particular QP charge-density profile) has been quite accurately calculated numerically hierarchy .

## III Haldane Hierarchy and Jain Sequence

Knowing that the QE–QE and QH–QH interactions are generally repulsive, Haldane proposed haldane  condensation of QE’s and QH’s into the hierarchy of “daughter” states at the series of Laughlin values of $`\nu _{\mathrm{QE}}`$ or $`\nu _{\mathrm{QH}}`$. In these states, the appropriate QP’s correlate with one another in the same way as the electrons do in the Laughlin states, and their elementary excitations are simply a new generation of QP’s. Assumming such QP condensation at each level of the hierarchy one would predict incompressibility of the whole electron system at all filling factors given by any odd-denominator fraction. This is in striking disagreement with the experiments, and the reason for this discrepancy is that although Coulombic, the QP–QP interaction at short range is not quite identical to the electron–electron repulsion responsible for the Laughlin correlations. As a result, the QP’s form Laughlin liquids only at very few of the Laughlin fractions, which eliminates all but a few valid “hierarchy” fractions hierarchy , in good agreement with the experiment. The same series of fractions arise naturally in the CF picture. These are the states at $`\nu =(2p+1/n)^1`$, corresponding to a number $`n`$ of completely filled CF LL’s each carrying flux $`2p\varphi _0`$.

The new FQH states pan  occur at the values of $`\nu `$ from outside the Jain sequence, and thus corresponding to only partial filling of a CF LL. Hence, in the CF picture, their incompressibility implies role of CF–CF interactions. In the QP hierarchy picture, the new states either coincide with the “invalid” fractions (e.g., $`\nu =4/11`$) or are new fractions altogether (e.g., $`\nu =3/8`$). In both cases it is clear that the origin of observed incompressibility lies in the special form of QP–QP correlations, and that these correlations are of a new (non-Laughlin) type.

## IV QP–QP Pseudopotential

The nature of QP correlations depends critically on the form of pseudopotential $`V_{\mathrm{QP}}()`$ describing their pair interaction energy $`V_{\mathrm{QP}}`$ as a function of relative pair angular momentum $``$. We have shown earlier parentage ; fivehalf  that the correlations are of the Laughlin type (i.e., the particles tend to avoid pair states with one or more of the smallest values of $`=1`$, 3, …) only if $`V()`$ is “superharmonic” at the relevant values of $``$ for a given filling factor $`\nu `$ (specifically, at $`=2p1`$ for $`\nu (2p+1)^1`$, where $`p=1`$, 2, …). Laughlin correlations defined in this way justify reapplication of the CF picture to the QP’s to select the lowest states of the whole many-body spectrum, and lead to the incompressible QP “daughter” states of the standard CF hierarchy sitko . The superharmonic repulsion is defined as one for which $`V`$ decreases more quickly than linearly as a function of the average particle–particle separation $`r^2`$ for the consecutive pair eigenstates labeled by $``$. In spherical geometry haldane , most convenient for finite-size calculations, this means that $`V`$ increases more quickly than linearly as a function of $`L(L+1)`$, i.e., of the squared total pair angular momentum $`L=2l`$, where $`l`$ is the single-particle angular momentum.

The qualitative behavior of the QP–QP interaction pseudopotential $`V_{\mathrm{QP}}()`$ at short range is well-known from numerical studies of small systems hierarchy ; sitko ; lee . On the other hand, the repulsive character of the QP–QP interaction and the long-range behavior of $`V_{\mathrm{QP}}()^{1/2}`$ follow from the fact that QP’s are charged particles (the form of QP charge density affects $`V_{\mathrm{QP}}`$ only at short range, comparable to the QP size). Combining the above arguments, it is clear that the dominant features of $`V_{\mathrm{QE}}`$ are the small value at $`=1`$ and a strong maximum at $`=3`$. Analogous analysis for the QH’s yields maxima at $`=1`$ and 5, and nearly vanishing $`V_{\mathrm{QH}}(3)`$.

## V QP Pairing

It is evident that $`V_{\mathrm{QE}}`$ does not support Laughlin QE–QE correlations. Instead, we expect that at least some of the QE’s will form pairs (QE2) at $`=1`$. A paired state would be characterized by a greatly reduced fractional parentage $`𝒢`$ parentage  from the strongly repulsive $`=3`$ state compared to the Laughlin correlated state, and have lower total interaction energy $`E=\frac{1}{2}N(N1)_{}𝒢()V()`$. Let us stress that such pairing is not a result of some attractive QE–QE interaction, but due to a tendency to avoid the most strongly repulsive $`=3`$ pair state. At sufficiently high QE density this can only be achieved by having significant $`𝒢(1)`$, which can be interpreted as pairing into the QE2 molecules. By analogy, the QH pairing is expected at $`=3`$. The range of $`\nu _{\mathrm{QP}}`$ at which pairing can be considered is limited by the condition that the separation between the pairs must exceed the pair size. While for the QE pairs this is satisfied at any $`\nu _{\mathrm{QE}}<1`$, the QH pairing can only occur at $`\nu _{\mathrm{QH}}<1/3`$.

## VI Laughlin Correlations Between Pairs

Having established that the QP fluid consists of (bosonic) QP2 molecules, the QP2–QP2 interactions need be studied to understand correlations. The QP2–QP2 interaction is described by an effective pseudopotential $`V_{\mathrm{QP}_2}()`$ that includes correlation effects caused by the fact that the two-pair wavefunction must be symmetric under exchange of the whole QP2 bosons and at the same time antisymmetric under exchange of any pair of the QP fermions. This problem is analogous to that of interaction between the electron pairs in the $`n=1`$ LL fivehalf .

Although we do not know $`V_{\mathrm{QP}_2}()`$ accurately, we expect that since it is due to the repulsion between the QP’s that belong to different QP2 pairs, it might be superharmonic at the range corresponding to the QP2–QP2 separation. Our preliminary numerical results for four QE’s seem to support this idea. However, in contrast to the $`n=1`$ electron LL fivehalf , the lack of accurate data for $`V_{\mathrm{QP}}`$ at the intermediate range makes such calculations uncertain.

## VII Condensed Pair States

The assumption of Laughlin correlations between the QP2 bosons implies the sequence of Laughlin condensed QP2 states that can be conveniently described using the “composite boson” (CB) model fivehalf . Let us use spherical geometry and consider the system of $`N_1`$ fermions (QP’s) each with (integral or half-integral) angular momentum $`l_1`$ (i.e., in a LL of degeneracy $`g_1=2l_1+1`$). Neglecting the finite-size corrections, this corresponds to the filling factor $`\nu _1=N_1/g_1`$. Let the fermions form $`N_2=\frac{1}{2}N_1`$ bosonic pairs each with angular momentum $`l_2=2l_1_1`$, where $`_1`$ is an odd integer. The filling factor for the system of pairs, defined as $`\nu _2=N_2/g_2`$ where $`g_2=2l_2+1`$, equals to $`\nu _2=\frac{1}{4}\nu _1`$. The allowed states of two bosonic pairs are labeled by total angular momentum $`L=2l_2_2`$, where $`_2`$ is an even integer. Of all even values of $`_2`$, the lowest few are not allowed because of the Pauli exclusion principle applied to the individual fermions. The condition that the two-fermion states with relative angular momentum smaller than $`_1`$ are forbidden is equivalent to the elimination of the states with $`_24_1`$ from the two-boson Hilbert space. Such a “hard core” can be accounted for by a CB transformation with $`4_1`$ flux quanta attached to each bosontheorem . This gives effective CB angular momentum $`l_2^{}=l_22_1(N_21)`$, LL degeneracy $`g_2^{}=g_24_1(N_21)`$, and filling factor $`\nu _2^{}=(\nu _2^14_1)^1`$.

The CB’s defined in this way condense into their only allowed $`l_2^{}=0`$ state ($`\nu _2^{}=\mathrm{}`$) when the corresponding fermion system has the maximum density at which pairing is still possible, $`\nu _1=_1^1`$. At lower filling factors, the CB LL is degenerate and the spectrum of all allowed states of the $`N_2`$ CB’s represents the spectrum of the corresponding paired fermion system. In particular, using the assumption of the superharmonic form of boson–boson repulsion, condensed CB states are expected at a series of Laughlin filling factors $`\nu _2^{}=(2q)^1`$. Here, $`2q`$ is an even integer corresponding to the number of additional magnetic flux quanta attached to each CB in a subsequent CB transformation, $`l_2^{}l_2^{}=l_2^{}q(N_21)`$, to describe Laughlin correlations between the original CB’s of angular momentum $`l_2^{}`$. From the relation between the fermion and CB filling factors, $`\nu _1^1=(4\nu _2^{})^1+_1`$, we find the following sequence of fractions corresponding to the Laughlin condensed pair states, $`\nu _1^1=q/2+_1`$. Finally, we set $`_1=1`$ for the QE’s and $`_1=3`$ for the QH’s, and use the hierarchy equation hierarchy , $`\nu ^1=2p+(1\pm \nu _{\mathrm{QP}})^1`$, to calculate the following sequences of electron filling factors, $`\nu `$, derived from the parent $`\nu =(2p+1)^1`$ state

$$\nu ^1=2p+1(2+q/2)^1,$$
(1)
where “$`+`$” corresponds to the QE’s and “$``$” to the QH’s. Remarkably, all the fractions reported by Pan et al. are among those predicted for the $`\nu =1/3`$ parent. Note also that the same values of $`q=1`$, 2, 4, and 8 describe both observed QE and QH states. This indicates similarity of the QE–QE and QH–QH pseudopotentials and suggests that both $`V_{\mathrm{QE}_2}`$ and $`V_{\mathrm{QH}_2}`$ may be superharmonic only at the corresponding four values of $``$ (in such case, the remaining fractions could not be observed even in most ideal samples).

## VIII Conclusion

We have studied the QP–QP interactions leading to novel spin-polarized FQH states in the lowest LL. Using the knowledge of QP–QP pseudopotentials and a general dependence of the form of correlations on the super- or subharmonic behavior of the pseudopotential, we have shown that QP’s form pairs over a certain range of filling factor $`\nu _{\mathrm{QP}}`$. Then, we argued that the correlations between the QP pairs should be of Laughlin type and proposed a hierarchy of condensed paired QP states. The proposed hierarchy of fractions agrees remarkably well with the recent experiment of Pan et al. pan .

###### Acknowledgements.

The authors thank Jennifer J. Quinn and Josef Tobiska for helpful discussions and acknowlege support from US DoE grant Grant DE-FG 02-97ER45657. AW acknowledges support from KBN grant 2P03B02424. KSY acknowledges support from KOSEF grant R14-2002-029-01002-0.


text:
# The Case for Asymmetric Dust Around a C-Rich AGB Star

## 1 Introduction

IRAS 15194–5115 (II Lup) was classified as a carbon star by Epchtein et al. (1987) and Meadows et al. (1987). At $`12\mu \mathrm{m}`$ it is the third brightest carbon star in the sky. It is a Mira variable for which Le Bertre (1992) gave an approximate period of 580 days. The circumstellar environment of the star has been studied by a number of workers. For instance Ryde et al. (1999) modelled the radio and far infrared rotational lines of CO, whilst Groenewegen et al. (1998) modelled the infrared spectral energy distribution. The star has similarities to the proto-type mass-losing carbon Mira IRC+10 216 (CW Leo) which has a pulsation period of 630 days. Ryde et al. derive a mass-loss rate of $`1\times 10^5\mathrm{M}_{}\mathrm{yr}^1`$ for II Lup and suggest that the wind characteristics of the object have not changed over the past few thousand years. However, they note that their data do not preclude variations in the mass-loss rate in the inner parts of the circumstellar material by a factor of three. Both Ryde et al. and Groenewegen et al. find that the dust-to-gas ratio in II Lup is about twice that in the envelope of CW Leo. However, the main difference between II Lup and CW Leo is in the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio. Ryde et al. estimate this as 5.5 for II Lup. Whereas Kahane et al. (1992) found $`44\pm 3`$ for CW Leo. One possibility, suggested by Ryde et al. for the relatively low ratio in II Lup, is that the star has recently gone through a period of hot-bottom burning.

Drops in the near infrared brightness of a carbon Mira on a time-scale long compared with the pulsation period, were first seen in the Galactic Mira R For (Feast et al. 1984). In that case a decline at visual wavelengths was also reported. These events were attributed to increased absorption in the line-of-sight. Further observations of R For were published by Le Bertre (1992). The most extensive investigation of this phenomenon is that of Whitelock et al. (1997) who discuss the near-infrared light curves of 11 large amplitude carbon variables from data covering time periods of from 9 to 22 years together with visual observations for some of them.

Here we present and discuss $`JHKL`$ (1.2, 1.6, 2.2 and 3.5 $`\mu `$m) observations of II Lup over a time period of about 18 years. The results are of particular interest for at least three reasons. (1) During this time period the object underwent an extended obscuration phase. (2) The pulsation variations are well covered, early and late in the time period involved. (3) At one point during the obscuration phase, infrared speckle observations (Lopez et al. 1993) provide estimates of the relative contributions of the starlight and the dust-shell emission in the infrared. These estimates are helpful in interpreting the data presented here.

## 2 Observations

The SAAO near infrared photometry of II Lup is listed in Table 1. The times of the observations are given as Julian Date (JD) from which 2440000 has been subtracted. The $`JHKL`$ measurements were all made with the MKII photometer on the 0.75m telescope at SAAO, Sutherland. They are on the current SAAO system as defined by the standards of Carter (1990). The photometry is accurate to better than $`\pm 0.03`$mag at $`JHK`$ and to better than $`\pm 0.05`$mag at $`L`$. Three of the early measurements were published by Lloyd Evans & Catchpole (1989). This is the star called WO48 in that paper, but note that two of the JDs were inadvertently transposed there.

## 3 Discussion

### 3.1 The Light Curves

Figure 1 shows the $`JHKL`$ light curves The figure also includes as crosses the measurements from Le Bertre (1992) as well as the observations from Epchtein et al. (1987), Meadows et al. (1987) and Groenewegen et al. (1993). There are likely to be differences between the various photometric systems. This is clearly seen at $`L`$ where an $`L^{}`$ filter was used by the ESO observers. There is no clear evidence of a significant difference between the SAAO $`J`$ data and that of the other observers. In any case the non-SAAO data are not used except to define the general shape of the light variations. The very bright $`L^{}`$ magnitude at JD 2448314 is from Groenewegen et al. (1993). Their $`JHK`$ measures do not appear unusual. They only report this single observation and we have not considered their $`L^{}`$ measurement in the later discussion.

In the following we concentrate on discussing the $`J`$ and $`L`$ light curves. Those at $`H`$ and $`K`$ are intermediate between these. Fig. 2 shows the fit to the $`J`$ data of a sinusoid with $`P_0=575`$ days and a long-term trend, modelled as sinusoid with $`P_1=12\times P_0`$; a small contribution is also made by the first harmonics of the two periods, $`P_0/2,P_1/2`$. The data stream is not long enough to determine the nature of the long-term trend and the fact that it can be represented in this way is probably fortuitous.

The 575 day period, which we take as the pulsation period, is well defined. This is only a slight revision of the period of 580 days given by Le Bertre (1992). Evidently this revised period satisfies the early observations (i.e. JD 2447000-2448000) when the object was faint at $`J`$ and also the later ones (JD 2451600-2452600) when it was bright at $`J`$. This apparent constancy of the period suggests that there has been no significant changes in the properties of the underlying star during the period of observation. The amplitudes of the pulsation are slightly smaller at the brighter epoch (by $`0.2`$ mag at $`J`$, $`0.4`$ mag at $`H`$, $`0.2`$ mag at $`K`$ and $`\stackrel{<}{}0.1`$mag at $`L`$). In view of the known irregularities of Mira light curves it is not clear whether these small changes have any important significance. However, they may be due to changes in the relative contributions of the starlight and the dust emission to the total flux at the early and late epochs (see below).

The near constancy of the mean $`L`$ brightness is striking compared with the large changes at $`J`$ (a range of more than two magnitudes). There is, however, a small increase in the brightness at $`L`$ (of about 0.1mag) between the early (faint) and later (bright) epochs.

It is evident from Fig. 1 that the mean brightness of II Lup at $`J`$ underwent an extended minimum lasting from at least JD $``$ 2445800 to JD $``$ 2451800. Despite the sparseness of the data at some epochs it is clear that the $`J`$ variations during this period were complex and that the curves shown in Fig. 2 provide only an approximate fit to the long term variations. There was evidently a steep decline in mean $`J`$ brightness near JD 2446000. Between JD $``$ 2447000 and JD $``$ 2449200 there appears to have been a further, slower, decline followed by a rise in brightness. A relatively steep decline followed at JD $``$ 2449600. From JD $``$ 2450500 to the end of the time period under discussion (JD 2452500), the object has been brightening at $`J`$.

It is of interest to compare the brightness variations of II Lup at $`J`$ with some theoretical models. There has been extensive work over a number of years on the connection between pulsation, mass-loss and dust formation. See, for instance, the summaries and references to earlier work by two particularly active groups in Sedlmayr & Winters (2000), Fleischer et al. (2000), Winters (2003), Höfner & Dorfi (1997), Sandin & Höfner (2003), and Andersen et al. (2003). This work has centred on the formation of shock waves and dust in the circumstellar environment. A paper by Winters et al. (1994) makes predictions about the long term light curves expected from models of this type. These predictions are sensitive to the C/O element abundance ratio adopted in the models, as was particularly emphasized by Fleischer et al. (2000). For relatively high values of the ratio ($``$ 1.8) no long term modulation of the light curve is predicted. However, for lower values (e.g. $`1.3`$) a long term periodicity is predicted, new dust being formed not every pulsation cycle but on a longer, periodic, time scale. Such a simple periodicity is not apparent in the behaviour of the $`J`$ brightness of II Lup just described. This result is similar to that obtained by Whitelock et al. (1997) for a number of other carbon Miras. These authors pointed out that the infrared light curves of R For were complex and did not conform to the model predictions. R Vol and R Lep were similar to R For.

There are some carbon variables, e.g. R Scl (Whitelock et al. 1997), which may show long-term periodic modulation in their infrared light curves, though this needs confirmation. Note that the carbon variable V Hya with a secondary (optical) period of $`6160`$ days (Mayall 1965, Knapp et al. 1999, Olivier et al. 2001) is thought to be in a binary orbit of this period. However, it would appear that in general carbon Miras do not show the long term periodic modulation predicted by some models. The long term variations of these stars are much less regular. It may be noted that the model of Winters et al. (1994) has been criticized by Groenewegen (1997) in its application to CW Leo because it requires a different distance for the object to fit the 11 $`\mu `$m visibility curve from that needed to fit the spectral energy distribution. We may conclude that whilst such theoretical models may well form the basis for further advances, they do not fit all the observations. The current models are based on spherically symmetrical systems. In the next section we give reasons for supposing that there were significant departures from spherical symmetry in the case of the II Lup during the time period under discussion.

### 3.2 Evidence for a Dust Cloud Model

In view of the apparent constancy of the pulsation period through the time period of the observations discussed above, we assume in the following that there has been no overall change in the underlying star. The tentative suggestion (Whitelock et al. 1997) that in the case of R For an apparent change in pulsation period might be associated with a change in mass-loss rate, needs further investigation.

Infrared speckle observations of II Lup were obtained by Lopez et al. (1993) on JD 2448434 when the object was in the obscuration minimum (see Fig. 1 or 2). These workers used their observations together with a spherically symmetrical model to estimate the relative contributions of the dust emission and starlight to the total flux at various wavelengths. In view of the departure from spherical symmetry which we are suggesting in the present paper, the numerical results derived by Lopez et al. have to be treated with caution. Nevertheless it seems safe to assume from their work that at this time almost all of the flux at $`L`$ came from the dust shell whilst almost all of the flux at $`J`$ was starlight. The near constancy of the mean $`L`$ flux throughout the period of our observations then suggests a constant dust shell. On the other hand the variations at $`J`$ suggest changes in the shell opacity in the line-of-sight. This evidently requires an asymmetric model.

To roughly quantify these results we adopt the proportions of flux from the star and the shell at different wavelengths as given by Lopez et al. Table 1 shows that (by chance) we have a $`JHKL`$ observation on the same day as their observation. The object was then around phase 0.16 in the pulsation cycle. Lopez et al. (see their fig. 6) interpret their results to show that at that time there was a negligible contribution of heated dust in the shell to the flux at $`J`$. Most of the $`J`$ flux was direct starlight though about 10 percent of the flux was scattered starlight. In contrast, at $`L`$ 80 percent of the flux was dust emission and only 20 percent starlight.

Towards the end of the time period under discussion, when the object was bright at $`J`$, it was observed at almost exactly the same pulsation phase as that of the speckle work (JD 2452457 see Table 1). Call this epoch B and the epoch of the speckle work A (these epochs are marked in Fig. 1). Suppose the variation in $`J`$ between A and B is due to a change in dust absorption. The speckle work indicates that at epoch A 15 percent of the flux at $`H`$ was from heated dust. Adopting this value and assuming that the flux from heated dust remained constant, it follows from the photometry at A and B that the effective absorption law of the dust is $`\mathrm{\Delta }J_{}/\mathrm{\Delta }H_{}`$ = 1.24, where the asterisk indicates the stellar component of the flux111Without allowing for the dust component at $`H`$ this ratio is 1.35. This is similar to the effective reddening law found for the carbon Mira R For (Feast et al. 1984) which gives 1.33 for this ratio. The uncertainty of this quantity is difficult to estimate, but may be about 0.1 in the case of both variables. This agreement with R For, which has much bluer $`JL`$ colours (see below), suggests that the conclusions of Lopez et al. (1993), assigning the bulk of the $`J`$ flux to the stellar contribution, are broadly correct. In connection with their discussion of circumstellar dust formation round carbon variables, Andersen et al. (2003) consider three types of carbon-rich particles using data taken from Rouleau & Martin (1991) and Jäger et al. (1998). The predicted values for the above ratio range from 1.36 to 1.85 for the three types of particle. The rough agreement with the values actually observed is probably as good as can be expected in view of both observational and theoretical uncertainties. The fact that the speckle work indicates that dust emission contributes a large fraction (50 percent) of the flux at $`K`$ at epoch A, makes it impractical to extend this type of calculation to wavelengths longer than $`H`$.

If 80 percent of the flux at $`L`$ is from dust emission at epoch A then $`L_{sh}=0.15`$, where the subscript “sh” denotes shell, and $`L_{}=+1.36`$. The work on R For (Feast et al. 1984) gave $`\mathrm{\Delta }J/\mathrm{\Delta }L4.0`$ as the effective reddening ratio of the circumstellar shell. This value is uncertain as it might be affected by (an unknown) shell-emission component at $`L`$. But note that the range covered by the infrared colours of R For used to determine the reddening law are bluer than those of II Lup (e.g. a range in $`(JL)_0`$ of 4.5 to 5.4 for R For compared with values for II Lup at epochs A and B of 8.0 and 6.0). Thus the contribution of dust emission to the $`L`$ flux is likely to be lower for R For than for II Lup. The three dust models considered by Andersen et al. (2003) predict the ratio $`\mathrm{\Delta }J/\mathrm{\Delta }L`$ to be in the range from 3.2 to 5.8. If the shell flux at $`L`$ has remained constant between epochs A and B we can use an adopted ratio $`\mathrm{\Delta }J_{}/\mathrm{\Delta }L_{}`$ to predict the brightening of the stellar component of $`L`$, and hence the expected total brightness at $`L`$. For adopted values of $`\mathrm{\Delta }J_{}/\mathrm{\Delta }L_{}`$ of 3, 4, 5 and 6 one then predicts the $`L`$ magnitudes at epoch B to be $`0.55,0.51,0.49`$ and $`0.47`$. The observed value is $`0.49`$. Considering the various uncertainties there is no disagreement with this range of models. We may conclude that the observations are consistent with the hypothesis of constant emission from the shell and a varying stellar component between epoch A and B due to variable dust absorption.

It may be noted that if we ascribed the increase in brightness at $`J`$ between epochs A and B to an increase in the emission flux from the dust, then the dust flux at $`J`$ would need to increase by about 6 magnitudes. A similar increase at longer wavelengths is then anticipated unless the temperature structure of the shell changed markedly. Any such change would have to take into account that the pulsation amplitudes are little affected by the rise in brightness. Although the $`J`$ flux at B would then be essentially all from the dust flux whilst at A it is essentially all starlight. In addition the total $`L`$ flux has to be kept nearly constant. We therefore discard this possibility.

We might anticipate that if the increase in brightness at $`J`$ between A and B were due to the dispersal of a spherical shell of dust, then the circumstellar dust emission would be reduced at the same time. An estimate of this reduction could be made if we knew the intrinsic colours of the underlying star. These are not known. However, the carbon Miras in the LMC discussed by Feast et al. (1989) are all relatively bright optically and $`JH`$ is likely to be only slightly affected by circumstellar reddening. These twenty Miras have $`\overline{JH}=1.23`$ with a range of individual values from 0.70 to 1.67. The five carbon Miras with periods greater than 350 days have $`\overline{JH}=1.39`$ and a range 1.21 to 1.67. The longest period Mira in this group has a period of 418 days and $`JH=1.30`$. The above values need to be corrected for the interstellar reddening of the LMC stars. This is expected to be small ($`E_{(JH)}0.03`$ from Feast et al. 1989). The optically bright (unobscured) carbon Miras in our Galaxy cover a similar range in $`(JH)_0`$ colours (Feast et al. 1982). One might expect the intrinsic colours of carbon Miras to increase with period. This is certainly the case for oxygen Miras. For the purposes of illustration we consider the consequences of adopting intrinsic colours $`(JH)_0`$ of 1.0 or 1.5 for II Lup. We also adopt an interstellar absorption of $`A_V=0.8`$ mag from Groenewegen et al. (1998). It is then straightforward to show that with 80 percent of the flux at $`L`$ from the shell at epoch A, there would be a decrease in the total flux at $`L`$ between A and B of 0.12 or 0.23 mag, assuming that the temperature of the dust emission remained constant. These figures account for both the expected decrease in the shell emission and the increase in starlight at $`L`$. In fact the decrease in obscuration at $`J`$ between A and B is presumably due to to the dust concerned moving away from the star. Thus the temperature of this dust will decrease between A and B and the actual decrease in brightness at $`L`$ will be greater than the figures just given. However, we in fact see a small increase ($`0.1`$ mag) in the total flux at $`L`$ between A and B which (as discussed above) can plausibly be attributed to the increase in direct starlight at the later epoch. These considerations suggest that we are not dealing with the dispersal of a complete dust shell. But with the clearing (due to expansion from the star) of a dust cloud of limited extent and in the line-of-sight. Note that for simplicity we refer to a dust cloud. However, in view of the relative complexity of the $`J`$ light curve, ejection in the line-of-sight was presumably taking place over a considerable length of time and at a variable rate.

### 3.3 General Discussion

The conclusion reached above, that the long-term modulation of the light curves of II Lup is due to the ejection of a dust cloud (or clouds) of limited extent in the line-of-sight, leaving the circumstellar dust emission virtually unchanged, is consistent with other data on carbon Miras. In the case of R For, it was argued by Whitelock et al. (1997) that an obscuration phase was due to a dust cloud which could not be a complete spherical shell. These authors also point out that the bolometric absolute magnitude of R Lep calculated from observations from the near- to the mid-infrared and using the Hipparcos parallax of the object, is fainter than anticipated theoretically. This conclusion is not changed if one adopts the revised parallax of R Lep recently published by Knapp et al. (2003). Whitelock et al. (1997) note that this problem can be overcome if the circumstellar emission from R Lep is from a non-uniform shell with higher than average absorption in the line-of-sight.

One, extreme, model for a non-spherical dust shell would be a circumstellar disc viewed edge-on. Disc models have been suggested to explain polarization observations of the carbon variable RW LMi (GL1403) (Cohen & Schmidt 1982) which has some similarities to R For (see Whitelock et al. 1997). Theoretical disc models have also been proposed (Dorfi & Höfner 1996). Such a model might be appropriate for the unusual object V Hya (e.g. Dorfi & Höfner 2000) which is probably a binary (see above). However, it is not certain that it is applicable to carbon Miras in general. For instance if the disc is optically thick to starlight we would expect significant changes in apparent bolometric luminosity between pole-on and edge-on views. For a disc with total opening angle (measured at the centre of the star) of 20 degrees this would amount to about two magnitudes. The scatter of long period carbon Miras with dust shells in the LMC about a bolometric period-luminosity relation is less than one magnitude (Whitelock et al. 2003) and much of this is likely to be due to the very poor temporal coverage of the stars at mid-infrared wavelengths.

In the sample of LMC carbon Miras studied by Whitelock et al. (2003) there are four which show evidence of obscuration phases. These, together with their $`K[12]`$ colours (outside major obscuration) are; IRAS 05300–6651, 7.29; TRM 72, 5.93; IRAS 05009–6616, 6.28; TRM 88, 4.44. The colours order the objects according to the relative contributions of the near- and mid-infrared to their bolometric fluxes. For the last three objects the near-infrared data are sufficient to allow a comparison between the bolometric magnitudes obtained using $`JHKL`$ values at the bright and faint phases of the obscurations. For these three stars (in the order given above) the differences between these values are only 0.1, 0.1 and 0.4 mag. The $`K[12]`$ colours of the last of these, TRM 88, are amongst the smallest of the long period sample of carbon Miras used by Whitelock et al. to define a period-luminosity relation and the object has a large obscuration range ($`\mathrm{\Delta }K`$ 1 mag). Thus the conclusion drawn above regarding the relative narrowness of the period-luminosity relation is not affected by the way the bolometric magnitudes were calculated by Whitelock et al.

The optically bright and relatively short period ($`P<420`$ days) carbon Miras in the LMC discussed by Feast et al. (1989) fit tight period-luminosity relations at $`K`$ and $`M_{bol}`$. Such stars thus appear to have thin shells and their $`K`$ and $`M_{bol}`$ magnitudes are not seriously affected by circumstellar reddening. However, the Galactic carbon Mira, R Lep, mentioned above is an interesting intermediate case. With $`K[12]=2.41`$ its near infrared flux will contribute more to a derived $`M_{bol}`$ than in the case of the obscured LMC Miras mentioned above. Thus higher than average shell absorption in the line-of-sight could well lead to a considerable underestimate of the bolometric luminosity as suggested by Whitelock et al. (1997). II Lup, itself, with $`K[12]=3.39`$ may also be in this category.

Further work on the obscured LMC stars might allow one to discriminate decisively between a model with a disc of large vertical extent and a quasi-spherical model. At present the latter seems the more promising for most of the stars. Both models require small scale variations due to dust clouds in order to explain object like II Lup unless one invokes a precessing disc model. Even then some non-uniform structure in the disc would be necessary to take into account the complexity of the variations at $`J`$.

There are a number of other observations which point towards non-uniform shell models. High resolution observations show that the inner region of the circumstellar material around CW Leo is highly structured and variable (e.g. Tuthill et al. 2000, Men’shchikov et al. 2001 and the summary in Monnier et al. 2003), whilst remarkable images of the outer region show complex arclet structures (de Laverny 2003). Furthermore CO observations show the shells of carbon variables (mostly classed as semi-regular) to be clumpy (e.g. Olofsson et al. 1996, 1998).

It is interesting to note that obscuration events appear to be relatively rare in single oxygen Miras (see Whitelock et al. 2000 and also Bedding et al. 2003). It is not clear whether or not this difference between carbon and oxygen Miras could be due to the type of particle involved. Whitelock (1987, 2003) has drawn attention to the fact that obscuration minima are relatively frequent in symbiotic systems containing oxygen Miras. She suggests that this might be due to the wind from the hot component in the system producing inhomogeneities in the distribution of the dust.

The long term modulation of the light curves of carbon Miras by dust clouds ejected in the line-of sight is similar in principle to the model for the obscuration minima of R Coronae Borealis (RCB) variables, with the star itself left unchanged (a good example of this in the case of the pulsating RCB star RY Sgr is given in Feast (1979) and reproduced in Clayton (1996)). If the $`\mathrm{C}_2`$ bands seen in emission during obscuration phases of some carbon Miras (Lloyd Evans 1997) come from the outer atmosphere of these stars they could be explained as becoming visible when the main body of the star was dimmed by a dust cloud. This is the model for the emission spectrum seen in RCB stars during minima (Feast 1979). Whitelock et al. (1997) note some similarities between the long term variations of R For and RCB stars. Some differences between the details of the Mira obscuration phases and those of typical RCB stars, are, however, to be expected in view of the major differences between the size and temperatures of the stars involved, the order of magnitude difference in the expansion velocities of the circumstellar material, and possible differences in the effective reddening law of the shell.

Understanding the formation of the dust is a formidable problem and it remains to be seen whether the theoretical advances cited above can be combined with an asymmetrical ejection (dust cloud) model. The basic cause of asymmetrical mass ejection may be non-uniformities in the stellar atmosphere. At least in the case of oxygen Miras, interferometry has indicated departures from circular symmetry which might be due areas of different brightness (and temperature) (e.g. Lattanzi et al. 1997), possibly large convection shells. It has also been suggested (Soker & Clayton 1999) that dust formation in both AGB and RCB stars takes place preferably above cool magnetic spots. Alternatively, instabilities in the dust forming region itself may lead to the formation of dust clouds (Woitke et al. 2000).

## 4 Conclusions

Superposed on the pulsational brightness variations of II Lup at $`J`$ there is a long term, large amplitude variation. This is complex and, as in the case of R For (Whitelock et al. 1997) and probably other stars, it is not easily fitted by available models which predict long-term periodic variations for mass-losing carbon Miras. These models are spherically symmetrical. Combining our $`JHKL`$ data with the high resolution work of Lopez et al. (1993) allows us to conclude that the effective reddening law of the circumstellar material at $`J`$ and $`H`$ is similar to that found for R For (Feast et al. 1984). The overall properties of the star itself and the radiation from the circumstellar dust have remained practically constant over the $`18`$ years of observation. The long-term variations at $`J`$ are attributed to a dust cloud or clouds ejected in the line-of sight. Taken together with other evidence the most likely model for objects of this kind is a quasi-spherical dust shell, formed in material moving away from the star, with small scale irregularities which can produce large variations in the line-of-sight absorption as they form and disperse. The model is then similar to that proposed for RCB variables which are also carbon-rich (though hydrogen deficient). It is hoped that future interferometric work will be able to test the model we propose for II Lup.

## Acknowledgements

This paper is based on observations made at the South African Astronomical Observatory (SAAO). We are grateful to the following who each made some of the observations: Brian Carter, Robin Catchpole, Dave Laney, Karen Pollard and Greg Roberts.


text:
# Coupled Ripplon-Plasmon Modes in a Multielectron Bubble

## I Introduction

Multielectron bubbles are 0.1 $`\mu `$m – 100 $`\mu `$m sized cavities inside liquid helium, that contain helium vapor at vapor pressure and a nanometer-thick electron layer anchored to the surface of the bubble VolodinJETP26 . The dependence of the radius of the bubble on the pressure applied to the bubble TemperePRL87  allows to vary the surface density of the spherical two-dimensional electron system in the bubble over nearly four orders of magnitude. As such, these objects offer the prospect to study the two-dimensional electron system on helium in regimes hitherto inaccessible. In particular, this would allow to study the two-dimensional electron Wigner crystal at densities not currently achievable for electrons on a flat helium surface.

Early experimental evidence for 2D Wigner crystallization of electrons on a liquid-He surface GrimesPRL42  relied on the detection of the coupled plasmon-ripplon modes FisherPRL42 . In Ref. TemperePRL87 , the ripplon modes characteristic for a bubble inside liquid helium were calculated and investigated as a function of the pressure applied on the bubble. In this communication we extend this investigation and take into account the redistribution of charge along the bubble surface when the bubble deforms. This allows us to derive the dispersion relation for the spherical ‘plasmon’ modes and the plasmon-ripplon coupling peculiar for the spherical surface.

## II Small amplitude deformations and charge redistribution

The surface of a bubble in liquid helium can be described by a function $`R(\mathrm{\Omega })`$ which gives the distance of the helium surface from the center of the bubble, in a direction given by the spherical angle $`\mathrm{\Omega }=\{\theta ,\varphi \}`$. In equilibrium (at zero pressure) the multielectron bubble has a spherical surface with radius $`R_\text{b}`$. We expand the deformation away from the spherical surface in spherical harmonics:

$$R(\mathrm{\Omega })=R_\text{b}+\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}Q_\mathrm{}mY_\mathrm{}m(\mathrm{\Omega }).$$
(1)
where $`R_\text{b}`$ is the angle-averaged radius of the bubble, $`Q_\mathrm{}m`$ is the amplitude of the spherical harmonic deformation mode indexed by $`\{\mathrm{},m\}`$, and $`Y_\mathrm{}m(\mathrm{\Omega })`$ is the corresponding spherical harmonic. These modes of deformation are referred to as ripplons in analogy with the surface modes of a flat surface of liquid helium. In Ref. TemperePRL87  the dispersion relation for the ripplons in a multielectron bubble has been studied as a function of pressure and related to the bubble stability. In this communication,we extend the results of TemperePRL87  in order to take the coupling and mixing of the ripplons and the electronic modes into account.

In a multielectron bubble, the electrons are confined to a thin (1 nm) layer at the inner surface of the bubble. This layer is anchored to the surface of the bubble, so that when the surface deforms, the layer conforms to the new bubble surface. The electrons can redistribute themselves inside the spherical layer, so that the surface density of electrons is no longer uniform. We describe the surface density with the function $`n_S(\mathrm{\Omega })`$ giving the number of electrons in a solid angle $`d\mathrm{\Omega }`$ around $`\mathrm{\Omega }:`$

$$n_S(\mathrm{\Omega })=\frac{N}{4\pi R_b^2}+\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}n_\mathrm{}mY_\mathrm{}m(\mathrm{\Omega }).$$
(2)
This function is normalized such that $`N=n_S(\mathrm{\Omega })R_b^2𝑑\mathrm{\Omega }`$ is the total number of electrons in the bubble or on the droplet. The $`n_\mathrm{}m`$’s represent the amplitudes of charge redistributions corresponding to spherical harmonics $`Y_\mathrm{}m`$.

The total potential energy of the multielectron bubble or droplet can be separated into several contributions: (i) a term from the surface tension energy $`U_S=\sigma S`$, with $`\sigma 3.6\times 10^4`$ J/m2 and $`S`$ the surface of the deformed bubble or droplet; (ii) a pressure related term $`U_p=pV`$, with $`p`$ the difference in pressure outside and inside the bubble and $`V`$ the volume of the bubble; (iii) a term representing the electrostatic energy $`U_C`$ of the electron layer. In this last term, the quantum corrections (such as the exchange contribution) to the electrostatic energy can be neglected ShungPRB45 . The first two terms, $`U_S`$ and $`U_p`$, have been derived in TemperePRL87 :

$`U_S`$ $`=`$ $`4\pi \sigma R_b^2+{\displaystyle \frac{\sigma }{2}}{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}\left(\mathrm{}^2+\mathrm{}+2\right)\left|Q_\mathrm{}m\right|^2,`$ (3)
$`U_p`$ $`=`$ $`{\displaystyle \frac{4\pi }{3}}pR_b^3+pR_b{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}\left|Q_\mathrm{}m\right|^2.`$ (4)
Whereas $`U_S`$ is an exact expression, $`U_p`$ is an expansion up to second order in $`Q_\mathrm{}m`$ and the third order term has been neglected. In what follows, we assume small amplitude deformations and small amplitude charge redistributions such that

$`\sqrt{\mathrm{}(\mathrm{}+1)}\left|Q_\mathrm{}m\right|`$ $``$ $`R_b,`$ (5)
$`\sqrt{\mathrm{}(\mathrm{}+1)}\left|n_\mathrm{}m\right|`$ $``$ $`{\displaystyle \frac{N}{4\pi R_b^2}}.`$ (6)
The electrostatic potential $`V(𝐫)`$ of the deformed MEB with a non-uniform surface electron density is calculated straightforwardly within the framework of the Maxwell equations by expanding the potential inside and outside the deformed bubble in their respective spherical decompositions and imposing the electrostatic boundary conditions at the surface. The potential energy associated with this electrostatic potential is given by

$$U_C=(1/2)n_S(\mathrm{\Omega })V(R,\mathrm{\Omega })𝑑\mathrm{\Omega }.$$
(7)
If we keep only terms up to second order in both $`n_\mathrm{}m`$ and $`Q_\mathrm{}m`$, we find

$`U_C`$ $`=`$ $`{\displaystyle \frac{e^2N^2}{2\epsilon R_b}}+2\pi e^2R_b^3{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left|n_\mathrm{}m\right|^2}{\epsilon _1\mathrm{}+\epsilon _2\left(\mathrm{}+1\right)}}`$ (8)
$`{\displaystyle \frac{e^2N^2}{8\pi \epsilon R_b^3}}{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\epsilon _1\mathrm{}^2\epsilon _2\left(\mathrm{}+1\right)}{\epsilon _1\mathrm{}+\epsilon _2\left(\mathrm{}+1\right)}}\left|Q_\mathrm{}m\right|^2`$
$`e^2N{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{}+1}{\epsilon _1\mathrm{}+\epsilon _2\left(\mathrm{}+1\right)}}n_\mathrm{}mQ_\mathrm{}m^{}.`$
where $`e`$ is the electron charge, $`\epsilon _1`$ is the dielectic constant of the medium inside the surface, and $`\epsilon _2`$ is the dielectric constant of the medium outside the surface. For a multielectron bubble, $`\epsilon _2=\epsilon `$ with $`\epsilon =1.0572`$ the dielectric constant of liquid helium, and $`\epsilon _1=1`$ (if the vapour pressure of the helium is low enough, the bubble has vacuum inside).

The ripplon contribution to the kinetic energy of the MEB is associated with the motion of the liquid helium surface. Following a derivation of Lord Rayleigh for oscillating liquid droplets RayleighPRS29 , we find

$$T_\text{r}=\frac{\rho R_b^3}{2}\dot{R}_b^2+\frac{\rho R_b^3}{2}\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{\mathrm{}+1}\left|\dot{Q}_\mathrm{}m\right|^2.$$
(9)
where $`\rho 145`$ kg/m3 is the density of liquid helium. Note that for a bubble, $`\mathrm{}+1`$ appears in the denominator instead of $`\mathrm{}`$ for the droplet in Lord Rayleigh’s treatment.

The ‘classical’ kinetic energy of the electrons is given by

$$T_\text{e}=\frac{m_\text{e}}{2}\underset{j=1}{\overset{N}{}}\dot{𝐰}^2(𝐫_j)$$
(10)
with $`m_\text{e}`$ the electron mass, and $`𝐰(𝐫_j)`$ represents the displacement of electron $`j`$ out of its equilibrium position $`𝐫_j`$. Since we assumed (6), we can use the formula

$$T_\text{e}=\frac{Nm_\text{e}}{8\pi }_{\text{surface}}\dot{𝐰}^2(𝐫)d^2𝐫$$
(11)
and express the field of displacements $`𝐰(𝐫)`$ as a sum of a longitudinal field $`𝐰_\text{L}(𝐫)`$ and transverse field $`𝐰_\text{T}(𝐫).`$ We investigate the effect of the longitudinal field, which can be written as a gradient of a scalar potential. The (divergence-free) transverse field is not considered here since it does not couple to the ripplons. Using

$$𝐰_𝐋(𝐫)=1\frac{4\pi R_b^2}{N}n_S(\theta ,\varphi )$$
(12)
we can express the kinetic energy of the electrons as

$$T_\text{e}=\frac{1}{2}\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\frac{4\pi m_\text{e}R_b^6}{N\mathrm{}(\mathrm{}+1)}\left|\dot{n}_\mathrm{}m\right|^2$$
(13)
We have checked that our approach yielded, for the 3D and the 2D electron gas, the known expressions for the plasma frequencies.

## III Results and discussion

The full Lagrangian of the bubble, including the ripplon modes and charge redistribution modes, is given by substituting expressions (9),(13),(3),(4),(8) in $`=T_\text{r}+T_\text{e}U_SU_pU_C`$. The result can be brought in the following form

$``$ $`=`$ $`_R+{\displaystyle \underset{\mathrm{},m}{}}{\displaystyle \frac{M_{\mathrm{}}}{2}}\left(\dot{Q}_\mathrm{}m^2\omega _r^2(\mathrm{})Q_\mathrm{}m^2\right)`$ (14)
$`+{\displaystyle \underset{\mathrm{},m}{}}{\displaystyle \frac{m_{\mathrm{}}}{2}}\left(\dot{n}_\mathrm{}m^2\omega _p^2(\mathrm{})n_\mathrm{}m^2\right)+{\displaystyle \underset{\mathrm{},m}{}}\gamma _{\mathrm{}}n_\mathrm{}mQ_\mathrm{}m.`$
This is the central result of this communication. We will now proceed to discuss the result term by term. The first term in expression (14) contains the Lagrangian describing the radial motion (the anharmonic breathing mode):

$$_R=\frac{\rho R_b^3}{2}\dot{R}_b^2\frac{4\pi }{3}pR_b^34\pi \sigma R_b^2\frac{N^2e^2}{2\epsilon R_b}.$$
(15)
In what follows we assume that $`\dot{R}_b=0`$ and the bubble radius is given by its equilibrium value.

### III.1 Ripplon modes

The Lagrangian $``$ also contains a part representing the harmonic oscillation of the ripplon modes, with oscillator mass term

$$M_{\mathrm{}}=\frac{\rho R_b^3}{\mathrm{}+1},$$
(16)
and bare ripplon frequency:

$$\omega _r(\mathrm{})=\sqrt{\frac{1}{M_{\mathrm{}}}\left[\sigma (\mathrm{}^2+\mathrm{}+1)+pR_b\frac{N^2e^2}{4\pi \epsilon R_b^3}\frac{\mathrm{}^2\epsilon (\mathrm{}+1)}{\mathrm{}+\epsilon (\mathrm{}+1)}\right]}$$
(17)
The bare ripplon frequencies and their dependence on the pressure were the subject of TemperePRL87 . Note that, in the limit of very large bubbles, $`R_b\mathrm{}`$, the dispersion relation for ripplons on a flat helium surface Atkins  is recovered: if the momentum $`k`$ is identified with $`\mathrm{}/R_b`$, then the typical $`k^{3/2}`$ is retrieved. Note furthermore that, if $`R_b`$ is set equal to the equilibrium coulomb radius, the modes $`\mathrm{}=1,2`$ obtain a zero frequency in agreement with the result of Salomaa and Williams SalomaaPRL47 . For typical multielectron bubbles with radii in the range of 1-100 $`\mu `$m (and numbers of electrons of the order $`10^410^7`$) the ripplon frequencies for $`\mathrm{}<1000`$ lie typically in the MHz-GHz range.

### III.2 Longitudinal plasmon modes

Next, the Lagrangean (14) contains a part representing the (small amplitude) harmonic oscillation of the classical charge redistribution around its equilibrium. The oscillator mass term for these oscillations is

$$m_{\mathrm{}}=\frac{4\pi m_\text{e}R_b^6}{N\mathrm{}(\mathrm{}+1)},$$
(18)
and the corresponding frequencies are given by

$$\omega _p(\mathrm{})=\sqrt{\frac{Ne^2}{m_\text{e}R_b^3}\frac{\mathrm{}(\mathrm{}+1)}{\mathrm{}+\epsilon (\mathrm{}+1)}}.$$
(19)
In the limit of large bubbles, this frequency corresponds to that of the longitudinal plasmon frequency of a 2D electron system on the helium surface PlatzmanPRB10 . If the electrons form a Wigner lattice, this is the longitudinal phonon frequency. If the momentum $`k`$ is identified with $`\mathrm{}/R_b`$, then the typical $`k^{1/2}`$ disperson for this longitudinal plasmon is retrieved. Based on this correspondence, we will refer to these modes in a multielectron bubble as the ‘(longitudinal) plasmon modes of the MEB’ (or ‘longitudinal phonon modes of the MEB’ if the electrons form a Wigner lattice).

For typical multielectron bubbles with radii in the range of 1-100 $`\mu `$m (and numbers of electrons of the order $`10^410^7`$) the longitudinal plasmon frequencies for $`\mathrm{}<1000`$ lie typically in the GHz-THz range. For bubbles with $`N>10^5`$, the longitudinal plasmon frequencies are larger than the ripplon frequencies, $`\omega _p(\mathrm{})\omega _r(\mathrm{})`$.

The equations of motion for the electron charge redistribution on a deformed bubble are given by:

$`{\displaystyle \frac{d}{dt}}{\displaystyle \frac{}{\dot{n}_\mathrm{}m}}`$ $`=`$ $`{\displaystyle \frac{}{n_\mathrm{}m}}`$ (20)
$``$ $`m_{\mathrm{}}\ddot{n}_\mathrm{}m=m_{\mathrm{}}\omega _p^2(\mathrm{})n_\mathrm{}m+\gamma _{\mathrm{}}Q_\mathrm{}m.`$
When $`\omega _p(\mathrm{})\omega _r(\mathrm{})`$ the electrons can redistribute much faster on the surface of the bubble than the bubble surface can deform. Thus, we can make an adiabatic approximation to find the charge redistribution of the electrons on a bubble with a given deformation. We find from (20) that if the bubble surface deformation is described by a given set of spherical deformation amplitudes $`Q_\mathrm{}m`$, the equilibrium charge distribution on the deformed bubble must satisfy the relation

$$\frac{n_\mathrm{}m}{n_0}=\frac{(\mathrm{}+1)}{2}\frac{Q_\mathrm{}m}{R_b},$$
(21)
with $`n_0=N/(4\pi R_b^2)`$.

### III.3 Ripplon - longitudinal plasmon coupling

The last term in the Lagrangean (14) represents the coupling between the longitudinal plasmons and the ripplons, with coupling strength given by

$$\gamma _{\mathrm{}}=Ne^2\frac{(\mathrm{}+1)}{\mathrm{}+\epsilon (\mathrm{}+1)}.$$
(22)
For electrons on a flat helium surface, such a coupling between ripplons and longitudinal plasmons/phonons was derived by Fisher et al. FisherPRL42  and detected experimentally GrimesPRL42 . For electrons on the inner surface of a deformed bubble, we find a similar coupling, but only ripplon and longitudinal plasmon modes which have the same angular momentum couple to each other. After the diagonalisation of the ripplon-plasmon part of the Lagrangian (14), we arrive at the eigenfrequencies,

$$\mathrm{\Omega }_{1,2}\left(\mathrm{}\right)=\sqrt{\frac{1}{2}\left[\omega _p^2(\mathrm{})+\omega _r^2(\mathrm{})\right]\pm \frac{1}{2}\sqrt{\left[\omega _p^2(\mathrm{})\omega _r^2(\mathrm{})\right]^2+4\gamma _{\mathrm{}}^2}}.$$
(23)
In Fig. 1, we show the eigenfrequencies $`\mathrm{\Omega }_1\left(\mathrm{}\right)`$ and $`\mathrm{\Omega }_2\left(\mathrm{}\right)`$ for $`N=10^5`$ as a function of $`\mathrm{}`$ and $`p.`$ The frequency $`\mathrm{\Omega }_2\left(\mathrm{}\right)`$ is close to the ripplon frequency derived for an MEB within the approximation of a conducting surface TemperePRL87 . Consequently, this frequency can be attributed to the ripplon modes, renormalized due to the ripplon-plasmon mixing. The admixture with the longitudinal plasmon mode is weak. The other branch of oscillations with the frequencies $`\mathrm{\Omega }_2\left(l\right)`$ can be related to the longitudinal plasmon mode admixed with a small component of ripplon nature. In typical multielectron bubbles, the difference in frequency between the ripplon and the plasmon modes weakens the coupling between these modes. Unlike for electrons on a flat surface FisherPRL42  the mixing is weak. Finally, note that the multielectron bubble will be stable when both $`\mathrm{\Omega }_1^2`$ and $`\mathrm{\Omega }_2^2`$ are positive. This condition is equivalent to

$$\omega _p^2(\mathrm{})\omega _r^2(\mathrm{})>\gamma _{\mathrm{}}$$
(24)
By substituting the results (19), (17) and (22) in the inequality (24) it can easily be seen that when the radius of the bubble is equal to the equilibrium radius, the ripplon-plasmon mixing does not change the criterion of stability for a MEB formulated in Ref. TemperePRL87 .

## IV Conclusions

In this communication we studied the longitudinal plasmon-ripplon modes in a multielectron bubble. The central result is the Lagrangean (14) with the ripplon and plasmon frequencies (17),(19) and the coupling strength (22). Unlike for electrons on a flat surface, we find that the mixing of the modes is weak for typical multielectron bubbles, because the bare longitudinal plasmon and ripplon frequencies are different ($`\omega _p(\mathrm{})\omega _r(\mathrm{})`$). The conditions for bubble stability discussed in Ref. TemperePRL87  remain unaltered by the ripplon-plasmon mixing. The analytic expressions for the longitudinal plasmon frequencies and the coupling strength as a function of the angular momentum allow us to derive a compact formula describing the charge redistribution taking place on a deformed bubble. For a bubble deformation in a particular spherical harmonic mode $`\{\mathrm{},m\}`$, the charge will redistribute itself according to the same spherical harmonic mode, and the amplitudes of deformation will be proportional to each other.

## Acknowledgements

Discussions with J. Huang are gratefully acknowledged. J. T. is supported financially by the FWO-Flanders. This research has been supported by the Department of Energy, Grant DE-FG02-ER45978, and by the GOA BOF UA 2000, IUAP, the FWO-V projects Nos. G.0071.98, G.0306.00, G.0274.01, WOG WO.025.99N (Belgium).