text:
# Interactions of vortices with rarefaction solitary waves in a Bose-Einstein condensate and their role in the decay of superfluid turbulence.

## I Introduction

Low temperature superfluids and recently discovered dilute Bose-Einstein condensates are often modelled by the Gross-Pitaevskii equation (GP), which is a nonlinear Schröedinger equation on the one-particle wave function, $`\psi `$. This equation provides a simple framework for study many fundamental hydrodynamical properties of condensates. In particular, vortex-sound interactions in condensates are receiving increasing attention over last couple of years vinen1 ; leadbeater1 ; br9 ; berloff ; vinen2 ; leadbeater2 . It has been suggested that the emission of sound by vortex reconnections and vortex motion is the only active dissipation mechanism responsible for the decay of superfluid turbulence. The decay of superfluid turbulence via Kelvin wave radiation and vortex reconnections was studied in the framework of the GP equation leadbeater2 . In this study the collision of two small vortex rings were analysed and the loss of the vortex line due to reconnections and Kelvin wave radiation was numerically evaluated.

The goal of this paper is to consider the effect of the vorticity-free solitary waves that together with vortices are created in Bose-Einstein condensates either during the self-evolution of a Bose gas from a strongly nonequilibrium initial state or evolved from local density depletions of a condensate. We show that the interactions with such waves trigger the loss of the vortex line via Kelvin wave radiation and that these interactions enhance the dissipation of the vortex tangle.

We write the GP equation in dimensionless form as

$$2\mathrm{i}\frac{\psi }{t}=^2\psi +(1|\psi |^2)\psi ,$$
(1)
in dimensionless variables such that the unit of length corresponds to the healing length $`\xi `$, the speed of sound is $`c=1/\sqrt{2}`$, and the density at infinity is $`\rho _{\mathrm{}}=|\psi _{\mathrm{}}|^2=1`$. The numerical integration was performed using a finite difference scheme. The faces of the computational box were open to allow sound waves to escape; this is achieved numerically by applying the Raymond-Kuo technique rk . In the turbulence simulations the periodic boundary conditions were used to conserve the total number of particles in the system.

In a seminal paper, Jones and Roberts jr  numerically integrated the GP equation (1) and determined the entire sequence of solitary wave solutions of the GP equation, such as vortex rings, vortex pairs, and finite amplitude sound waves named rarefaction pulses. They showed the location of the sequence on the momentum, $`p`$, energy, $``$, plane. In three dimensions they found two branches meeting at a cusp where $`p`$ and $``$ assume their minimum values, $`p_m`$ and $`_m`$. As $`p\mathrm{}`$ on each branch, $`\mathrm{}`$. On the lower branch the solutions are asymptotic to large vortex rings.

As $``$ and $`p`$ decrease from infinity along the lower branch, the solutions begin to lose their similarity to large vortex rings. Eventually, for a momentum $`p_0`$ slightly greater than $`p_m`$, they lose their vorticity ($`\psi `$ loses its zero), and thereafter the solitary solutions may better be described as ‘rarefaction waves’. The upper branch consists entirely of these and, as $`p\mathrm{}`$ on this branch, the solutions asymptotically approach the rational soliton solution of the Kadomtsev-Petviashvili Type I equation.

The Jones-Roberts (JR) solitons are the only known disturbances that propagate with a constant velocity. Notice, however, that there many other waves that change their velocity and shape during their motion. For instance, a strong perturbation of a rarefaction pulse on the upper branch of the dispersion curve causes it to collapse onto the lower branch and become a vortex ring. During its transition the wave looses its energy and momentum and the minimum of its density decreases gradually before it reaches zero. In the next stage of the evolution the wave becomes a non-axisymmetric ring which radius increases until it reaches an axisymmetric solitary state on the lower branch of JR dispersion curve. All these intermediate time snapshots of the wavefunction of the condensate have higher energy than any solitary wave on the lower branch below the final axisymmetric vortex ring state although they may have the same minima of the density in the transition. In what follows we reserve the term “rarefaction pulse” or “rarefaction wave” to describe a finite amplitude sound wave which is a JR solitary wave. The scenario just described illustrates a typical mechanism in which the vortex line length may be created and increased as the result of perturbations of rarefaction pulses. In pade  we have also considered a creation of vortex rings as a result of energy and momentum transfer between two interacting rarefaction pulses. But given a complex tangle of interacting vortices we expect that the loss of the vortex line length will dominate its creation to account for the experimentally observed decay of superfluid turbulence.

In pade  we have developed an algorithm for finding approximations to the JR solitary wave solutions with the correct asymptotic behaviour at infinity. An axisymmetric solitary wave moving with the velocity $`U`$ along the $`x`$axis is accurately approximated by $`\psi (x,s)=1+u(x,s)+\mathrm{i}v(x,s)`$ where

$`u`$ $`=`$ $`{\displaystyle \frac{a_{00}+a_{10}x^2+a_{01}s^2+mc_{20}^{7/4}U(2x^2f(U)s^2)}{(1+c_{10}x^2+c_{01}s^2+c_{20}(x^2+f(U)s^2)^2)^{7/4}}},`$
$`v`$ $`=`$ $`x{\displaystyle \frac{b_{00}+b_{10}x^2+b_{01}y^2mc_{20}^{7/4}(x^2+f(U)s^2)^2}{(1+c_{10}x^2+c_{01}s^2+c_{20}(x^2+f(U)s^2)^2)^{7/4}}},`$ (2)
where $`a_{ij},b_{ij},c_{ij}`$ and the dipole moment $`m`$ can be determined from the series expansion and are functions of $`U`$. Also in (2), $`s^2=y^2+z^2`$, $`x=xUt`$, and $`f(U)=12U^2`$. Notice, that (2) represents a vortex ring as well as a vorticity-free rarefaction pulse depending on $`U`$. If $`2<a_{00}<1`$, then (2) represents a vortex ring as the power series expansion around zero shows; if $`1<a_{00}<0`$, then (2) gives an approximation of a rarefaction pulse, and $`a_{00}=1`$ is a borderline case, such that the solitary wave has a single zero of the wavefunction and, therefore, can be called a point defect. These approximations can be used as initial conditions in the numerical simulation that study interactions, when the initial state is prepared by multiplying the wavefunctions of the distant individual solitary waves. Without an accurate starting point in numerical calculations it would be impossible to separate clearly the effect of interactions from the evolution of each solitary wave by itself as it settled down from a poor initial guess.

## II Creation of rarefaction pulses in condensates

Here we study the excitations created by the collisions between rarefaction pulses and vortices and the effect these collision have on the evolution of vortex line length in a regime of superfluid turbulence with a tangle of vortices. We start by describing three typical scenario in which a large amount of rarefaction pulses is created, so that their collisions with vortices do become significant. The process of strongly nonequilibrium BEC formation in a macroscopically large uniform weakly interacting Bose gas was elucidated in svist  using numerical integration of the GP equation. As the system evolves from weakly turbulent state to state of strong turbulence, the phases of the complex amplitudes of the wave field $`\psi `$ become strongly correlated and the period of their oscillations become comparable with the evolution times of the occupation numbers. This signifies that the quasicondensate is formed with appearance of a well-defined tangle of quantized vortices and localized vorticity-free solitary structures such as rarefaction pulses manifesting the start of the final stage in the Bose gas evolution: the decay of superfluid turbulence. The right panel of Fig. 5 presented in svist  shows a single vortex ring surrounded by many rarefaction pulses as a result of the turbulent decay of the initial vortex tangle, which also implies that many more rarefaction pulses were present at the beginning of the decay.

Next we describe how the JR solitons and in particular rarefaction pulses can be created in Bose-Einstein condensates. It is generally believed that to create vortices it is necessary to transmit angular momentum by rotationally stirring the condensate with a laser beam. In bubble  we demonstrated that the collapse of a stationary spherically symmetric bubble can lead to the vortex nucleation. After the condensate fills the cavity, it begins to expand with growing density oscillations. These ’dips’ of the density are themselves unstable and the development of this instability leads to the creation of localized disturbances: vortex rings and rarefaction pulses.

Similarly, the evolution of a depletion of the condensate by laser beam after the laser was turned off may lead to creation of JR solitons. Notice, that there is no need to add angular momentum or deplete the density of the condensate to zero. We demonstrate this by considering a depletion of the condensate, such that the wavefunction of the condensate just before the laser is switched off is given by

$$\psi (𝐱)=\frac{1}{2}+\frac{1}{2}\mathrm{tanh}\left[0.01(x^2+0.1(y^2+z^2)100)\right].$$
(3)
This density depletion roughly corresponds to the experiments in the sodium condensates, see for instance raman , with the healing length $`\xi =0.3\mu m`$ and the Gaussian beam waist of more than $`10\xi `$. The oblate spheroidal form of the depletion (3) would be formed if such a laser beam was guided in the condensate in a short straight streak of the length of about $`50\xi `$.

We integrated the GP equation (1) using (3) as our initial field in the box of the volume $`V=120^3`$. The faces of the computational box were open to allow sound waves to escape. Fig. 1 shows the time snapshots of the density cross-sections for $`z=0`$. Twelve vortex rings of various radii and four rarefaction pulses moving outward from the center of the depletion are clearly seen at $`t=73`$ snapshot. Fig.2 gives the isoplots of the density at $`|\psi |^2=0.3`$. Rarefaction pulses are seen as small oblate spheroids.

The nucleation of vortices in a uniform condensate has been linked to critical velocities of the flow frisch ; br7 ; br8 . It has also been pointed out rica  that a moving object of a size comparable with the healing length generates rarefaction pulses rather than vortex rings when the velocity on the surface of this object exceeds the local speed of sound. Therefore, it is possible to generate the rarefaction pulses by guiding small objects with supercritical velocities through the condensate.

## III Collisions with a straight line vortex

In our first example we consider a collision between a rarefaction pulse and the straight-line vortex of unit winding number. The rarefaction pulse is moving with a constant self-induced velocity $`U=0.63`$ and belongs to the lower branch of the JR dispersion curve. The sequence of density isosurface plots illustrating the collision for two impact angles ($`\pi /2`$ and $`\pi /4`$) is shown in Fig. 3. Initially the center of the rarefaction pulse is 20 (Fig. 3(a)) or 10 (Fig. 3(b)) healing lengths apart from the axis of the vortex. During the collision the rarefaction pulse creates a distortion on the vortex line by exciting two Kelvin wave packets. In the close vicinity of the vortex line, the rarefaction pulse acquired vorticity it previously did not have and became a vortex ring (Fig. 3(a) at $`t=25`$ and $`t=30`$ and Fig. 3(b) at $`t=20`$), so we would expect that the result of the collision is similar to what was found in leadbeater2  for collisions of large vortex rings with small rings.

The dispersion relationship for Kelvin wave in the GP model in the limit $`k0`$ is roberts02

$$\omega \frac{1}{2}k^2(\mathrm{ln}(k)+0.003187\mathrm{}),$$
(4)
which gives a group velocity of the Kelvin wave packet as

$$U_g=\frac{\pi }{\lambda }\left(2\mathrm{ln}[\frac{\lambda }{2\pi }]1.00637\right),$$
(5)
where $`\lambda `$ is the dominating wavelength. Our calculations show that $`\lambda 8`$ healing lengths and $`U_g0.2`$ ($`U_g0.3c`$), where $`c`$ is the speed of sound) which fits (5) quite well. This result also agrees (when different scaling of the GP equation is taken into account) with the wavelength of the excited Kelvin wave during the vortex reconnections in leadbeater2 .

The collision with a rarefaction pulse from the upper branch of the JR dispersion curve excites a Kelvin wave packet of a similar central wavelength. Fig. 4 shows the contour plots of the density cross-sections of two such collisions for two different rarefaction pulses that belong to lower (left panel) and upper (right panel) branches of the JR dispersion curve. The wavelength of the created excitation is best seen by the intersections of zeros of real and imaginary parts of the wavefunction $`\psi `$ shown by solid and dashed lines correspondingly. Similarly, a collision with an offset generates a Kelvin wave packet of approximately the same dominating wavelength but a smaller amplitude.

We can estimate the amplitude of the Kelvin wave generated by a distant  rarefaction wave by reducing the problem to two dimensions. We shall assume that initially the straight line vortex is positioned at $`(0,0,z)`$ and the rarefaction wave is moving in the positive $`x`$direction along $`y=d_y<0`$ with $`\xi d_y`$. Under these assumptions, the wavelength of the perturbation along the vortex line is much larger than its amplitude, $`\eta `$, (the maximum displacement of the vortex center from the origin in $`xy`$plane), therefore, we can assume that the $`x`$ and $`y`$ components of the force exerted by a rarefaction wave on the vortex dominate the $`z`$component (in other words, the vortex may be considered as rigid). Also, if the incident angle of the rarefaction pulse and the direction of the vortex axis is nonzero, the amplitude calculated below would have to be multiplied by the cosine of this angle. In $`xy`$plane the vortex is advected by the rarefaction pulse, so that it instantaneous velocity coincides, to the leading order, with the velocity of the rarefaction pulse at the vortex location. The center of the rarefaction pulse is moving with the self-induced velocity, $`𝐔=(U,0,0)`$ and is advected by the vortex. The motion of the vortex and of the rarefaction pulse is, therefore, given by the system of the coupled ordinary differential equations

$`{\displaystyle \frac{dx^r}{dt}}`$ $`=`$ $`U{\displaystyle \frac{y^ry^v}{(x^rx^v)^2+(y^ry^v)^2}},`$
$`{\displaystyle \frac{dy^r}{dt}}`$ $`=`$ $`{\displaystyle \frac{x^rx^v}{(x^rx^v)^2+(y^ry^v)^2}},`$ (6)
$`{\displaystyle \frac{dx^v}{dt}}`$ $`=`$ $`V_x(x^vx^r,y^vy^r),`$
$`{\displaystyle \frac{dy^v}{dt}}`$ $`=`$ $`V_y(x^vx^r,y^vy^r),`$
where $`(x^r(t),y^r(t))`$ and $`(x^v(t),y^v(t))`$ are the positions of the centres of the rarefaction pulse and the vortex, correspondingly, at time $`t`$. $`V_x`$ and $`V_y`$ are the velocity components induced by the rarefaction pulse on the vortex given by $`V_x=S/x`$ and $`V_y=S/y`$, where $`S`$ is the phase of the wavefunction $`\psi =R\mathrm{exp}^{iS}`$, so that

$$dS=\frac{(1+u)udvvdu}{(1+u)^2+v^2},$$
(7)
where the real, $`1+u`$, and imaginary, $`v`$, parts of $`\psi `$ are approximated by (2). The initial conditions are given by $`(x^v(0),y^v(0))=(0,0)`$ and $`(x^r(0),y^r(0))=𝐝=(d_x,d_y)`$. Fig. 5 shows the trajectories of the vortex for $`U=0.63`$ and $`𝐝=(10,10),(10,20),(20,10)`$ and $`(20,20)`$. It is clear from the figure, that the maximal amplitude of the displacement, $`\eta `$, does not depend strongly on $`d_x`$, but on $`d_y`$. Fig. 6 plots $`\eta `$ as a function of $`d_y`$ for two rarefaction pulses that belong to the lower ($`U=0.63`$) and upper ($`U=0.69`$) branches of the JR dispersion curve. For large offsets, $`d_y`$, ($`d_y>30`$) the maximum amplitude, $`\eta `$, decays exponentially with $`d_y`$, so that $`\eta d_y^{1.24}`$ for $`U=0.63`$ and $`\eta d_y^{1.62}`$ for $`U=0.69`$.

## IV Collisions with a vortex ring

Next we consider the collision of a rarefaction pulse with a vortex ring. As in the previous cases, such a collision excites Kelvin wave with release of sound energy. Fig. 7 illustrates the time snapshots of the collision between the vortex ring of radius $`R=5.1`$ moving along the $`x`$axis and the rarefaction pulse moving along the line $`y=R`$. The total vortex line is dramatically increased at time $`t11`$ (from the line length $`\mathrm{}32`$ to $`\mathrm{}50`$) due to the rarefaction pulse being transformed into a small vortex ring during a short time interval, $`\delta t5`$. After that the ring loses its energy to Kelvin wave radiation and stabilises to axisymmetric vortex ring of slightly smaller radius at $`t100`$, so that the total loss of the vortex line is about $`3\%`$. The length of the vortex line does not change for $`t>100`$. We repeated the calculations for the rarefaction pulse from the upper branch of the JR dispersion curve and found that the ring stabilises after losing $`4\%`$ of its length. These results may seem to contradict the findings of leadbeater2 , where the length of the vortex line continues to decay after the collision. The calculations of leadbeater2  were done in computational box with periodic boundaries, so the sound waves generated during the reconnections of a large ring with a small ring continue to excite Kelvin waves on the vortex line with a continuous loss of the vortex line to sound radiation. The rate of the decay of the vortex line in this case should depend on the size of the computational box. In our calculations the faces of the computational box were open to allow sound waves to escape rk , which allowed us to compute the effect of a single collision with the finite amplitude sound wave.

For very small vortex rings ($`R1`$) the collision with a rarefaction pulse causes all the energy and momentum to be converted into sound waves with a complete loss of the vortex line.

## V Decay of a vortex tangle

Next, we estimate the effect that the interactions with rarefaction pulses have on a vortex tangle. This time we perform the simulations similar to those of kivotides ; leadbeater2 , in which the vortex tangle was created by four colliding rings in a periodic box. But in our simulations we regulate the initial amount of sound waves in our computational box by introducing them as randomly oriented rarefaction waves in addition to the vortex rings.

For simplicity we shall assume that all sound waves initially consist of three types of rarefaction pulses moving in random directions. The approximations to these solitary waves, that were developed in pade , become very useful in the construction of our initial conditions. If the rarefaction pulse given by (2) is rotated by angle $`\alpha `$ to the $`x`$axis in the $`(x,y)`$ plane and by angle $`\beta `$ to the $`y`$axis in the $`(y,z)`$ plane, then its wavefunction is given by $`\psi (x^{},s^{})`$ where $`x^{}=x\mathrm{cos}\alpha y\mathrm{sin}\alpha `$, $`s^{}=\sqrt{(y^{}\mathrm{cos}\beta z\mathrm{sin}\beta )^2+(z\mathrm{cos}\beta +y^{}\mathrm{sin}\beta )^2}`$, and $`y^{}=y\mathrm{cos}\alpha +x\mathrm{sin}\alpha `$. The angles of rotation $`\alpha `$ and $`\beta `$ as well as the center of the rarefaction pulse are chosen randomly.

Our calculations were performed in a periodic box with volume, $`V=80^3`$, which is divided by $`161^3`$ grid points with a spacing $`0.5`$. The numerical scheme is 4th order globally accurate with 4th order Runge-Kutta time integration with the time step $`0.025`$. A grid spacing of $`0.25`$ together with the time step $`0.00625`$ was also used to test the accuracy of the numerical method. Initially there are four vortex rings of the radius $`R=30.1`$ centered at $`(R,1,0),(R,1,0),(1,R,0)`$ and $`(1,R,0)`$ that are moving towards the center of the box. $`N_{rare}`$ rarefaction pulses are distributed randomly in the internal half of the box away from the vortices. The results for the dynamics of the vortex line length for three sets of computations with $`N_{rare}=0,100`$ and $`200`$ are plotted in Fig. 8. The density isosurfaces plots for $`N_{rare}=100`$ are given in Fig. 9. In the case of nonzero number of rarefaction pulses there is a rapid initial growth of the length of the vortex line due to the evolution of some of the rarefaction pulses into small vortex rings as the result of the energy transfer among them pade . Consequently, the length of the vortex line shows the balance between vortices being created and vortex line being destroyed due to reconnections and Kelvin wave radiation. Our simulations show that apart from the initial growth of the vortex line in the simulations with $`N_{rare}=100`$ initially, the dynamics for $`N_{rare}=100`$ and $`N_{rare}=0`$ is very similar for in the time interval $`[400,2600]`$. On the other hand, the case with $`N_{rare}=200`$ is quite different. During the time interval $`[500,1300]`$, the decay of the vortex tangle is much faster than for two other runs. Whereas after $`t1800`$ it is the run with $`N_{rare}=200`$ that shows a slower decay. We speculate that this can be best explained by the following. During the moderate times, the vortex line loss due to the excitations of the vortex line created by rarefaction pulses depends on how many rarefaction pulses survived as their energy is gradually being converted to high-frequency waves. The more of rarefaction pulses were present in the initial state, the longer at least some of them are present in the system, therefore, the state with $`N_{rare}=200`$ shows a faster decay during the moderate times. At larger times, all of the rarefaction pulses are destroyed by the collisions and vortices live on the background of high-frequency waves that limit their mobility and reduce the number of reconnections which slows down the decay of the vortex line length, with the total energy carried by these high-frequency waves being higher for the simulations with $`N_{rare}=200`$ initially.

For even larger times the situation is reversed one more time as reconnections are now infrequent and the main decay mechanism is through Kelvin wave radiation. Therefore, the rate of decay of vortex line is propotional to the amount of high-frequency waves (slowest decay for $`N_{rare}=0`$, fastest decay for $`N_{rare}=200`$.)

For Kelvin-wave cascade, where energy is transfered to a much shorter wavelengths with a cutoff below a critical wavelength, the vortex line density $`L=\mathrm{}/V`$ can be described by the Vinen equation vinen

$$\frac{dL}{dt}=\frac{\kappa }{2\pi }\chi _2L^2,$$
(8)
where $`\kappa =2\pi `$ in our dimensionless units and $`\chi _2`$ is a dimensionless coefficient. More accurately $`\chi _2`$ is a weak (logarithmic) function of $`L`$ and of other parameters such as the Kelvin cutoff and temperature. The logarithmic dependence on $`L`$ can be easily obtained in the context of the local induction approximation. The presence of a large amount of rarefaction pulses changes this weak dependence on $`L`$ and the decay of the vortex line approaches an exponential decay instead:

$$\frac{dL}{dt}=\frac{\kappa }{2\pi }\chi _2L.$$
(9)
In Fig. 10 we plotted the curve (the solution of (9))

$$\mathrm{}=\mathrm{}_0\mathrm{exp}(\chi _2t)$$
(10)
for $`\chi _2=0.00035`$ and $`\mathrm{}_0=760`$ to illustrate the exponential decay of the vortex line length.

Finally, we note that the presence of rarefaction waves in a condensate with regular vortex lattices has destabilising effect and facilitates reconnections. If vortices are antiparallel, the Kelvin wave created on the filaments as the result of interaction with rarefaction waves facilitates the growth of the Crow instability br9  that leads to vortex reconnections that destroy a significant amount of vortex line (as we estimated in br9 , the minimum fractional line loss of the pair of antiparallel vortices that are initially the distance $`h`$ apart and perturbed by a wave of the wavenumber $`k1`$ is about $`1\sqrt{kh/2}`$). The interactions of the rarefaction waves with the lattice of parallel vortices (created, for instance, in rotating condensates) also destabilises the array and creates a tangle as Fig. 11 illustrates.

## VI Conclusions

In summary, we established several mechanisms for creation of the rarefaction pulses in condensates. We have studied the effects of their interactions with straight line vortices, vortex rings and tangles of vortices. We showed that there are two competing mechanisms of these interactions. Firstly, the interactions of waves at close distances lead to the increase of the vortex line as the rarefaction pulses in the regions of lower density (in the vicinity of the vortex core or each other) may acquire circulation and become vortex rings. Secondly, rarefaction waves excite Kelvin waves on the vortex filaments causing a loss of the vortex line due to sound emission. At high vortex line densities our simulations suggest that the Kelvin wave radiation is enhanced by the presence of rarefaction pulses and this accounts for a dramatic increase in the rate of the decay of the vortex line length.

###### Acknowledgements.

The author acknowledges the support from the NSF grant DMS-0104288. The author is very greateful to Professor Paul Roberts for many discussions about the structure of rarefaction pulses and the importance of sound emission during her tenure at UCLA. The author is indebted to Professor Boris Svistunov for many useful comments about this manuscript.


format: text

text:
# An observation of a nascent fractal pattern in MD simulation for a fragmentation of an fcc lattice

## Abstract

To seek for a possible origin of fractal pattern in nature, we perform a molecular dynamics simulation for a fragmentation of an infinite fcc lattice. The fragmentation is induced by the initial condition of the model that the lattice particles have the Hubble-type radial expansion velocities. As time proceeds, the average density decreases and density fluctuation develops. By using the box counting method, it is found that the frequency-size plot of the density follows instantaneously a universal power-law for each Hubble constant up to the size of a cross-over. This cross-over size corresponds to the maximum size of fluctuation and is found to obey a dynamical scaling law as a function of time. This instantaneous generation of a nascent fractal is purely of dynamical origin and it shows us a new formation mechanism of a fractal patterns different from the traditional criticality concept.

thanks: Present address: Department of Theoretical Studies, Institute for Molecular Science, Myodaiji, Okazaki 444-8585, Japan

Since the pioneering work of Mandelbrot  a large effort has been directed towards understanding the ubiquitous manifestation of fractal structure in nature. Since fractal structure is intimately related to a scaling property, a customary way of identifying the fractal is to check if the power-law relation of the frequency as functions of the size holds. Malcai et al.  in this way have listed the fractal dimensions and the corresponding scaling ranges for a variety of physical systems. They pointed out that the finiteness of the fractal range is a typical feature for physical systems which is in sharp contrast with an ideal mathematical concept of fractal.

Let us therefore settle our target as “why and how a power-law relation in frequency-size relation appears frequently in nature no matter what its range”. It seems natural to adopt some fractal growth models to challenge this subject or to adopt some kinds of cellular-automata approaches. In these models, the power-law in frequency-size relation is established at an asymptotic time. Now the question is if we can detect a bud of fractal pattern in the short-time regime long before the asymptotic time. Anatomy of this formation process at short time is the subject of this paper.

In a typical fractal growth model of diffusion limited aggregation (DLA) , the short-time behavior is not studied directly. Its behavior is studied in the simulation study  and experimental study  of kinetic roughing model, in the model for fracture process  and in several cellular-automata approaches used to study the self-organized criticality . These models in short-time regime are, however, concentrated on the phenomena of dynamical scaling relation  where the time-dependence of some growth parameter is investigated. The behavior of the growth parameter, however, is not an sufficient answer to the growth mechanism of fractal pattern asymptotically achieved.

In order to investigate the short-time behavior of the formation of fractal pattern, we propose a new simulation model for the fragmentation of an infinite fcc lattice. For this purpose, we adopt a box-counting method  to search for the fragmentation process instead of the fragment mass distribution. It is because the mass distribution is only defined for the low-density limit where a fragment is separated each other, which is not the case for the short-time regime where the density is not low enough to identify the fragments. In fact, we will see later that the box counting method is an appropriate method to investigate the bud of fractal, which starts as far as the system starts to expand from its saturation density.

The dynamics of the fragmentation is caused in our model by the Hubble-like radial expansion that is imposed as the initial condition of the molecular dynamics simulation . The Hubble expansion is the simplest and most symmetric driving force for the fragmentation since any point in the system is equivalent to any other point, a perfectly local equivalence for everywhere when we deal with an infinite system. Infinite system is chosen to eliminate the effect of surface and thus we concentrate on a fragmentation initiated by the fracture occurring in the volume. We perform a molecular dynamics simulation using the Lennard-Jones (L-J) potential with a cut-off length $`r_c=2.5`$ . As the initial condition, the particles are located at the lattice points of the fcc lattice in the ground state, however, 10% of randomly selected lattice points are removed so that the initial density $`\rho `$ becomes 0.95 in a reduced unit (all quantities in what follow are measured in the reduced unit). This imperfection is introduced to prepare enough number of different initial states of the fragmentation. Our main results are essentially independent of the vacancy fraction, in the range 5% to 15%. The initial velocity of a particle at $`𝐫`$ is given by $`h𝐫`$ where $`h`$ is a dimensionless parameter. Thus the system mimics a uniformly expanding Hubble universe. When solving the equation of motion for the infinite system represented by one primitive cell and 26 replica cells, we impose a generalized periodic boundary condition , under which a replica cell with the origin at $`𝐑_i`$ moves with the velocity $`h𝐑_i`$. This means that the time evolution is completely determined by the one parameter $`h`$.

A complete nascent fractal structures emerge only for suitably selected value of $`h`$. When $`h`$ is too small, the effect of the attractive force dominates the expansion and one single droplet whose mass equals to the total number of particles in a unit cell is formed. Since this structure depends on the unit cell size, we judge it unphysical. When $`h`$ is too large, the expansion dominates the response of the system and the solutions to the equation of motion result in a simple enlarged copy of the initial configuration. Since the initial configuration with near saturation density is almost uniform, the resultant structure has trivial dimension 3, even after the expansion. It is only with $`h`$ in an intermediate range, namely $`0.1h0.4`$ that our numerical simulations result in an emergence of intriguing fractal structures.

In the box-counting method, a unit cubic cell of the size $`R_{\mathrm{cell}}`$ of expanding matter is divided into smaller cubic boxes of the size $`\delta =R_{\mathrm{cell}}/n`$ where $`n`$ is a positive integer. We count the number $`N(d)`$ of small boxes that contain at least one particle. The fractal dimension $`D_f`$ (box counting dimension) is introduced with the relation $`N(\delta )=R_{\mathrm{cell}}^3/\delta ^{D_f}`$, which is rewritten $`\mathrm{log}[N(\delta )]=3\mathrm{log}[R_{\mathrm{cell}}]D_f\mathrm{log}[\delta ]`$. When $`\mathrm{log}[N(\delta )]`$ is proportional to $`\mathrm{log}[\delta ]`$, $`D_f`$ becomes constant and the relation between $`N(\delta )`$ and $`\delta `$ satisfies a power-law. In the following discussion, we will use the $`\mathrm{log}\mathrm{log}`$ plot of $`N(\delta )`$ and $`\delta `$ to obtain $`D_f`$.

Figure 1 shows the time evolution of snapshots of the fractal structures and the results of the box counting method applied to these structures. The left column corresponds to $`h=0.1`$ and the right column to $`h=0.2`$. For each column, the time evolution proceeds from the top to the bottom. As time proceeds the density decreases. The three graphs presented are for the average densities $`\rho =`$0.20, 0.10, and 0.05. For our studies, we prepare 50 initial configurations each of which is composed of 1,235 particles (fcc lattice with 10% vacancies). $`N(\delta )`$ is calculated as an average of these 50 configurations. Note that for $`\delta <1`$, $`N(\delta )`$ is always equal to the total number of particles because two particles cannot come closer than $`\delta 1`$ due to the repulsive term of the L-J potential. From lines fitted by the method of least squares and from the location of the intersection point between the two lines we conclude that the physical size of the fractal structure grows with decreasing density and that $`D_f`$ does not change significantly during the time evolution for a given $`h`$, that is, $`1.95D_f1.98`$ for $`h=0.1`$ and $`1.66D_f1.74`$ for $`h=0.2`$. The value of $`D_f`$ for $`h=0.1`$ is larger than that for $`h=0.2`$ and the range of scale where $`D_f`$ remains constant for $`h=0.1`$ is wider than for $`h=0.2`$. We checked if these results are independent of the number of particles in the unit cell and found that a cell composed of 9,878 particles leads to almost the same results. In the snapshot picture, we show an example of a configuration with 9,878 particles.

A similar behavior is observed in the case of more rapid expansion. In Fig.2, we show snapshots of the fractal evolution and the corresponding plot of the results of the box counting method for $`h=0.3`$ and $`h=0.4`$ at $`\rho =0.05`$. From Figs.1 and 2, we conclude that a larger $`h`$ gives a smaller $`D_f`$ and narrower range of scale. Note that the intersection point of the 3-dimensional and $`D_f`$ lines corresponds to the maximum size of the density fluctuation in the system. Beyond this crossover point, the distribution of matter is uniform and its dimensionality is three.

A summary of fractal structure for $`h=0.1`$ and $`h=0.2`$ is shown in Fig.3 where the fitted lines used to determine $`D_f`$ in Fig.1 are depicted together with additional lines for $`\rho =0.40`$. For each figure, the four lines corresponding to $`\rho =0.40,0.20,0.10,0.05`$ are almost overlapping for small values of $`\delta `$. The intersection points of these lines with the four non-overlapping 3-dimensional lines in Fig.3a and Fig3b are shifted toward larger $`\delta `$ values with decreasing density. The overlap of the four lines for small value of $`\delta `$ is the most important result of this study. It means that there is a nontrivial physical process that tunes the box counting so regularly. This feature is in sharp contrast with the results for $`h\mathrm{}`$ as shown in fig.3c. In this case, the effect of expansion dominates the system response and the density snapshot is a trivial enlarged copy of the initial configuration that simply results in a parallel shift of the 3-dimensional lines in the figure. It is remarkable that the fractal structures born already at densities near $`\rho =0.4`$ are maintained at lower densities. The approximate constancy of $`D_f`$ and the overlap of $`D_f`$ lines as functions of time clearly shows that there is a special correlations governing the whole system. Although we have not yet fully understood the mechanism behind this phenomenon, it is clearly caused by the critical balance between the externally imposed expanding motion and the internal potential.

When we follow the time evolution beyond $`\rho =0.05`$, the least square fit of two lines to the data becomes worse. This degradation occurs when the internal expansion of connected aggregations of particles almost stops. After this time, only the space between the fragments expands. This is another reason that we used the word critical balance since we need an on-going competition between the expanding external field and internal system forces in order to get persistent fractal behavior. As a result, the range of scale for which we get a good least-square fit for the fractal dimension $`D_f`$ is limited to about one-order of magnitude.

Finally, we will see the meaning of the cross-over point. As was discussed before, the box-size of this crossover point is a measure of the maximum size of the density fluctuation at a given time. The time-dependence of this size is shown in Fig.4, where the cross-over box-size $`\delta _{cross}`$ for $`h=0.1`$ and for $`h=0.2`$ are plotted together with the cell size $`R_{\mathrm{cell}}`$ as a function of the quantity $`1+ht`$. The calculated points are shown by various symbols corresponding to fours system densities, $`\rho =`$0.4,0.2,0.1 and 0.05 from the left to right, and lines are fitting to these four points. The linear dependence of $`R_{\mathrm{cell}}`$ on the abscissa is the original assumption of boundary motion. On the other hand, the fitting to cross-over size yields two curves

$`\delta _{\mathrm{cross}}`$ $`=1.222(1+ht)^{2.6883}`$ $`\mathrm{for}h=0.1,`$ (1)
$`\delta _{\mathrm{cross}}`$ $`=1.252(1+ht)^{2.0436}`$ $`\mathrm{for}h=0.2.`$ (2)
These power-law dependence of the quantity on time is known as the dynamical scaling phenomena . @ The dynamical exponent depends sensitively on the expansion speed and is larger for slower expansion. Their nature seems to be quite different from the corresponding quantities examined so far . We like to stress here that when we extrapolate to t=0, the value of $`\delta _{cross}`$ becomes about 1.2, which nearly coincides with the smallest size of the system fluctuation. We can say, in other word, that the maximum size of the fluctuation starts to increase in coincidence with the start of expansion. The bud of a fractal starts to grow from the saturation density already.

In conclusion, we have shown that a set of many-body coupled Newton equations together with an externally imposed initial condition of uniform expansion for an fcc lattice directly leads to the birth and evolution of fractal structures for specific ranges of the initial expansion rate. The fractal dimension obtained with the box counting method remains constant up to a cross-over scale , which obeys the dynamical scaling law as a function of time. The most important finding is that the nascent fractal is generated from the very beginning of the expansion, which might give a hint to the original motivation of this study, “why the fractal-like structure prevails in nature?” .

We thank S. Chiba for encouragement and support and P. Möller for critical readings of the original manuscript. S.C. is grateful to supports from a Research fellowship of the Japan Society for the Promotion of Science for Young Scientists.


text:
# Energy Crisis in the Superbubble DEM L 192 (N 51D)

## 1 Introduction

Massive stars inject large amounts of mechanical energy into the ambient interstellar medium (ISM) via fast stellar winds and supernova explosions. The stellar energy feedback alters the physical conditions of the ISM and consequently may regulate the formation of future-generation stars. To gain a better understanding of the physical structure and evolution of the ISM in a galaxy, it is necessary to investigate empirically the conversion efficiency of stellar mechanical energies into interstellar kinetic and thermal energies.

Massive stars usually form in groups such as OB associations, and their fast stellar winds and supernova ejecta collectively sweep up the ambient ISM to form superbubbles (Mac Low & McCray, 1988). Superbubbles provide excellent laboratories to study the effects of stellar energy feedback on the ISM because the kinetic and thermal energies in a superbubble are well-defined and can be measured with little ambiguity. Superbubbles in the Large Magellanic Cloud (LMC) at 50 kpc are particularly well suited for such studies, as their stellar content and ISM can be observed in great detail with little foreground extinction.

We have chosen the superbubble DEM L 192 (Davies et al., 1976) in the LMC for a detailed study of the interaction between its stars and the ambient ISM. DEM L 192, also known as N 51D (Henize, 1956), hosts two OB associations, LH 51 and LH 54 (Lucke & Hodge, 1970). The stellar content of these two OB associations and the dynamics of the superbubble have been analyzed in detail by Oey (1996a) and Oey & Smedley (1998). They conclude that LH 51 and LH 54 are both young, $`3\times 10^6`$ yr old. The diffuse X-ray emission from the hot interior of DEM L 192 was detected by Einstein and ROSAT observations (Chu & Mac Low, 1990; Dunne et al., 2001), but high-quality spectra of the X-ray emission have been provided only by XMM-Newton observations. The emission line features in the X-ray spectra suggest an overabundance in oxygen and neon which is used to argue for the occurrence of supernova explosions (Bomans et al., 2002).

This paper reports our investigation of the energetics of the superbubble DEM L 192. In §2 we analyze archival XMM-Newton observations of DEM L 192 to determine the physical properties of the hot gas in the superbubble interior. In §3 we determine the thermal energy of the hot superbubble interior and the kinetic energies of the ionized and neutral gas in the superbubble shell. In §4 we compute the stellar wind energy and supernova explosion energy expected from the known stellar content of LH 51 and LH 54. Finally, a discussion of the overall energy budget of the DEM L 192 superbubble and our conclusions are given in §5.

## 2 XMM-Newton Observations of Hot Gas in DEM L 192

DEM L 192 was observed with the XMM-Newton Observatory on 2001 October 31 and November 1 using the EPIC/MOS1, EPIC/MOS2, and EPIC/pn CCD cameras (Obs. ID: 7194). The two EPIC/MOS cameras were operated in the Prime Full-Window Mode for a useful exposure time of 31.6 ks, and the EPIC/pn camera was operated in the Extended Prime Full-Window Mode for a useful exposure time of 24.6 ks. The thin filter was used for all observations.

We retrieved the XMM-Newton pipeline products from the XMM-Newton Science Archive (XSA)111 The XMM-Newton Science Archive is supported by ESA and can be accessed at http://xmm.vilspa.esa.es/xsa. and further processed the data using the XMM-Newton Science Analysis Software (SAS ver. 5.3.3) and the calibration files from the Calibration Access Layer available on 2002 December 13. The event files were screened to eliminate bad events (e.g., events due to charged particles) so that only events with CCD patterns 0-12 (similar to ASCA grades 0-4) were selected for the EPIC/MOS observations, and events with CCD pattern 0 (single pixel events) were selected for the EPIC/pn observation. We also screened the event files to remove periods of high background. To assess the background rate, we binned the counts over 50 s time intervals for each instrument in the 10-12 keV energy range which is dominated by the background. Only for a short period of time did the high background, with count rates $``$0.3 counts s-1 for the EPIC/MOS or $``$0.6 counts s-1 for the EPIC/pn, affect the observations; events from this high-background period were discarded. The resulting net exposure times are 30.6 ks and 23.2 ks for the EPIC/MOS and EPIC/pn observations, respectively.

### 2.1 Spatial Analysis

We extracted images from the EPIC/pn and EPIC/MOS observations in the 0.25–2.5 keV energy band, which includes most of the diffuse X-ray emission from DEM L 192. These images were then combined to obtain a higher S/N ratio image and to reduce the null exposure in the gaps between CCDs (the satellite was not dithered during the observations). The resulting raw EPIC image is presented in Fig. 1a. This image was further adaptively smoothed using Gaussian profiles with FWHM ranging from 2″ to 15″ and divided by a normalized exposure map to remove the vignetting of the telescope and sensitivity variations of the instruments. A grey-scale presentation of this adaptively smoothed EPIC image is shown in Fig. 1b. For comparisons, an H$`\alpha `$ image of the same region is displayed in Fig. 1c, and the X-ray contours extracted from the smoothed EPIC image are overplotted on the H$`\alpha `$ image in Fig. 1d.

Diffuse X-ray emission from the interior of DEM L 192 is clearly detected. The brightest X-ray emission region is coincident with the OB association LH 54, peaking particularly at the bright Wolf-Rayet star HD 36402. The X-ray peaks coincident with bright stars may be dominated by stellar emission; therefore, emission from these peaks is excluded from our analysis of diffuse X-rays. The diffuse X-ray emission shows a rough limb-brightened morphology in the main body of the superbubble, with the brightest diffuse X-ray emission projected interior to the brightest H$`\alpha `$ emission regions of the shell. Faint diffuse X-ray emission is detected to the north and southwest of the main body; however, these faint X-ray extensions do not necessarily imply blowouts, as they are bounded by faint H$`\alpha `$ filaments externally. The superbubble DEM L 192 must be in an inhomogeneous ISM, and the variations in the H$`\alpha `$ surface brightness reflect the density variations in the ambient ISM. The rough correspondence between the diffuse X-ray surface brightness and the nebular H$`\alpha `$ surface brightness suggests that a mixing of the dense, cooler, superbubble shell gas into the hot interior has occurred.

### 2.2 Spectral Analysis

The X-ray spectrum of DEM L 192 is extracted from only the event file of the EPIC/pn camera because its sensitivity is greater than that of the EPIC/MOS cameras. The source aperture, marked in Fig. 1a, includes a large elliptical region that covers the main body of DEM L 192 and the north extension and a smaller square region that covers the southwest extension; these two regions together encompass all diffuse X-ray emission within the superbubble. The background contribution was assessed using an aperture exterior to DEM L 192 with similar area to the source aperture. All point sources were excised from the source and background apertures. A circular region that includes most stars in LH 54 is also excluded from the source aperture to minimize the stellar contribution to diffuse emission. After the background subtraction, a total of $``$16,000 counts are detected in the source aperture. Some residuals of the background subtraction are noticeable at the energy of the Al K$`\alpha `$ fluoresce line at $``$1.5 keV. The energy range 1.45-1.55 keV, thus, has been omitted in the subsequent spectral analysis.

The background-subtracted EPIC/pn spectrum of DEM L 192, shown in Fig. 2, peaks at 0.5–0.6 keV, with a bright shoulder at $``$0.7 keV and a weaker peak at 0.9–1.0 keV. Above 1.0 keV, the spectrum is featureless with emission detected up to 3.5 keV. The overall spectral shape indicates thermal plasma emission, with spectral features corresponding to the He-like triplet of O vii at $``$0.57 keV, the H-like O viii line at $``$0.65 keV, and the He-like triplet of Ne ix at $``$0.92 keV. Therefore, we have adopted the MEKAL optically-thin plasma emission model (Kaastra & Mewe, 1993; Liedahl, Osterheld, & Goldstein, 1995) with 1/3 $`Z_{}`$ abundances as appropriate for the LMC and the absorption cross-sections from Balucińska-Church & McCammon (1992). The spectral fits are carried out using XSPEC ver. 11.2.0 by folding the model spectrum through the EPIC/pn response matrices and comparing the modeled spectrum to the observed spectrum in the 0.3–4.5 keV energy range. The $`\chi ^2`$ statistics are used to determine the best-fit model.

The best-fit MEKAL model with 1/3 $`Z_{}`$ abundances is very unsatisfactory (reduced $`\chi ^2`$ = 7.2), failing to fit the spectrum at 1–3 keV and showing large residuals at the Ne ix lines at $``$0.92 keV. Adding a second thermal component at a higher temperature can improve the spectral fits significantly (reduced $`\chi ^2`$ = 2.2), but the temperature of the second component is pegged at the upper limit allowed by XSPEC, i.e., $`kT`$ = 80 keV. A thermal component at such a high temperature is not expected in the ISM; furthermore, it is dominated by bremsstrahlung emission of which the spectral shape is featureless and similar to that of a power law. Therefore, we have added a power-law component, instead of a second thermal component, and the reduced $`\chi ^2`$ is lowered to 2.2 for the best-fit model. The large residuals at the Ne ix lines remain, and they can be minimized only if the Ne abundance is allowed to vary. By increasing the Ne abundance by a factor of 3 (i.e., solar Ne abundance and 1/3 solar for the other elements), we are able to obtain a best-fit model with a reduced $`\chi ^2`$ of 1.4, as shown in Fig. 2. This confirms the enhanced Ne abundance suggested by Bomans et al. (2002); however, we do not find compelling evidence for an enhanced O abundance, which has also been suggested by Bomans et al. (2002).

In our best-fit model for the diffuse X-rays from DEM L 192, the thermal component dominates the emission below 1 keV with a temperature $`kT`$ of 0.21$`\pm `$0.01 keV, an absorption column density of (3.4$`\pm `$1.7)$`\times `$1020 H-atom cm-2, and a normalization factor222The normalization factor $`A`$ is defined to be $`A`$=$`1\times 10^{14}N_\mathrm{e}^2dV/4\pi d^2`$, where $`d`$ is the distance, $`N_\mathrm{e}`$ is the electron number density, and $`V`$ is the volume. $`A=(1.6\pm 0.4)\times 10^3`$ cm-5. The power law component dominates the emission from 1 to 3 keV and has a photon index $`\mathrm{\Gamma }`$ of $`1.3\pm 0.2`$. The observed flux of diffuse X-ray emission from DEM L 192 is $`1.1\times 10^{12}`$ ergs cm-2 s-1 and the X-ray luminosity is $`4\times 10^{35}`$ ergs s-1 in the 0.3–3.0 keV band. The thermal emission accounts for $``$70% of the observed luminosity (but $``$85% of the total number of detected photons).

## 3 Energy Content of DEM L 192

For a superbubble like DEM L 192, we expect the stellar mechanical energy injected into the ISM to be stored in the thermal energy of its hot interior and the kinetic energy of the dense, swept-up shell. As the ionization front might be trapped in the swept-up shell, both ionized and neutral gas in the superbubble shell must be considered. Below we examine separately the thermal energy of the hot interior and the kinetic energies of the ionized shell gas and associated H i gas.

### 3.1 Thermal Energy of the Superbubble Interior

The thermal energy of the hot gas in a superbubble is $`E_{\mathrm{th}}=(3/2)kTNϵ_1V`$, where $`T`$ is the plasma temperature, $`N`$ is the total particle number density, $`ϵ_1`$ is the volume filling factor of the hot gas, and $`V`$ is the volume interior to the superbubble. Assuming a canonical helium to hydrogen number density ratio of $`N_{\mathrm{He}}/N_\mathrm{H}`$ = 0.1, the total particle number density, approximated by the sum of electron density, hydrogen density, and helium density, is $`NN_\mathrm{e}+N_\mathrm{H}+N_{\mathrm{He}}1.92N_\mathrm{e}`$. For a given volume filling factor $`ϵ_1`$, the rms $`N_\mathrm{e}`$ can be derived from the normalization factor $`A`$ of the best-fit model: $`N_\mathrm{e}`$ = $`2\times 10^7d(\pi A/ϵ_1V)^{1/2}`$. For a distance of $`d`$ = 50 kpc, a normalization factor of $`A=(1.4\pm 0.4)\times 10^3`$ cm-5, and a volume $`V=4.4\times 10^{61}`$ cm3 (approximated by a $`120\times 120\times 200`$ pc3 ellipsoid), we find $`N_\mathrm{e}=(0.031\pm 0.004)ϵ_1^{1/2}`$ cm-3. The total mass of the hot gas in the superbubble is $`1.3\times 10^3ϵ_1^{1/2}`$ M, and the total thermal energy of the hot gas is $`E_{\mathrm{th}}=(1.3\pm 0.3)\times 10^{51}ϵ_1^{1/2}`$ ergs. The filling factor $`ϵ_1`$ is most likely between 0.5 and 1.0. For $`ϵ_1=0.5`$, $`N_\mathrm{e}`$ = $`0.044\pm 0.004`$ cm-3, the hot gas mass is 920 M, and $`E_{\mathrm{th}}=(9\pm 2)\times 10^{50}`$ ergs.

The cooling timescale of the hot gas $`t_{\mathrm{cool}}`$ is $`(3/2)NkT/\mathrm{\Lambda }(T)`$, where $`\mathrm{\Lambda }(T)`$ is the cooling function for gas at temperature $`T`$ in units of ergs cm-3 s-1. For a plasma temperature of a few $`\times `$ 106 K, $`\mathrm{\Lambda }(T)/N_\mathrm{H}^2=7\times 10^{24}`$ ergs cm3 s-1 for 1/3 Z abundances (Dalgarno & McCray, 1972). With $`N_\mathrm{H}0.83N_\mathrm{e}`$, $`N1.92N_\mathrm{e}`$, and $`kT=0.21\pm 0.01`$ keV, the cooling timescale $`t_{\mathrm{cool}}`$ is $`200`$ Myr. Since the cooling timescale is much larger than the age of the OB associations LH 51 and LH 54, $``$ 3 Myr, only a negligible amount of thermal energy has been radiated away.

### 3.2 Kinetic Energy of the Ionized Superbubble Shell

The mass of the ionized superbubble shell can be derived from the H$`\alpha `$ luminosity, $`L_{\mathrm{H}\alpha }`$. The H$`\alpha `$ flux from DEM L 192 is $`f_{\mathrm{H}\alpha }=(860\pm 80)\times 10^{12}`$ ergs s-1 cm-2 (Kennicutt & Hodge, 1986). For a distance of $`d=50`$ kpc, $`L_{\mathrm{H}\alpha }=4\pi d^2f_{\mathrm{H}\alpha }=(2.6\pm 0.3)\times 10^{38}`$ ergs s-1. In addition, for ionized gas at temperatures of $``$104 K, $`L_{\mathrm{H}\alpha }=3.56\times 10^{25}N_\mathrm{e}N_\mathrm{p}Vϵ_2`$ ergs s-1, where $`N_\mathrm{e}`$ and $`N_\mathrm{p}`$ are the particle number densities of electrons and protons in units of cm-3, respectively, $`V`$ is the volume of the region in units of cm3, and $`ϵ_2`$ is the volume filling factor. Again, we approximate the superbubble volume by a $`120\times 120\times 200`$ pc3 ellipsoid, and the volume is $`V=4.4\times 10^{61}`$ cm3. Helium is most likely singly ionized, thus $`N_\mathrm{e}=1.1N_\mathrm{p}=(4.5\pm 0.5)ϵ_2^{1/2}`$ cm-3. The total mass of the superbubble shell is $`M=1.27N_\mathrm{e}Vϵ_2m_\mathrm{H}`$, where $`m_\mathrm{H}`$ is the mass of a hydrogen atom. Using the rms $`N_\mathrm{e}`$ derived from $`L_{\mathrm{H}\alpha }`$, we find the superbubble shell mass to be $`M=(4.2\pm 0.5)\times 10^{38}ϵ_2^{1/2}`$ g = $`(2.1\pm 0.2)\times 10^5ϵ_2^{1/2}`$ M. The expansion velocity of the DEM L 192 superbubble has been measured in the H$`\alpha `$ line to be $`v_{\mathrm{exp}}35`$ km s-1 (Meaburn & Terrett, 1980). Thus, the kinetic energy in the ionized superbubble shell is $`E_{\mathrm{kin}}=\frac{1}{2}Mv_{\mathrm{exp}}^2`$ = $`(2.6\pm 0.3)\times 10^{51}ϵ_2^{1/2}`$ ergs. If the thickness of the ionized superbubble shell is 1/10 of the radius, then the filling factor of shell gas is $`ϵ_2=0.27`$, the ionized shell mass is $`M1\times 10^5`$ M, and the ionized shell kinetic energy is $`E_{\mathrm{kin}}1\times 10^{51}`$ ergs.

### 3.3 Kinetic Energy of Associated H i Gas

The kinetic energy of the neutral gas in the swept-up shell can be derived from observations of the 21-cm line emission from the H I gas. Such observations of DEM L 192 have been extracted from H I surveys of the LMC made with the Australia Telescope Compact Array (ATCA) and the Parkes Observatory 64-m telescope, where single-dish observations from the Parkes telescope have been combined with the ATCA data to “fill-in” low spatial frequencies not observed by the interferometer. The ATCA and Parkes observations are detailed in Kim et al. (1998) and Kim et al. (2003), respectively.

The kinematics of the H I gas can be determined from the position-velocity plots in Fig. 3. The H I gas shows a broad velocity profile centered at the systemic velocity of $`V_{\mathrm{hel}}300`$ km s-1 and with a width of $``$70 km s-1, corresponding to twice the expansion velocity measured in H$`\alpha `$ emission, but no clear expanding shell structure can be identified. The H I channel maps (Dunne et al., 2004) and interstellar absorption observations of stars in DEM L 192 (Wakker et al., 2004) indicate that the H I velocity structure is rather complex and some velocity components may not be physically associated with the superbubble.

A total 21-cm line flux density of 470 Jy is measured from the combined radio observations within the dotted polygon outlined in the H I column density map in Fig. 3. Given the LMC distance of 50 kpc, we estimate a total H I mass of 3 $`\times `$ 105 M from this flux density. The total kinetic energy of the H I gas can be approximated by 3/2 of the second moment of the integrated velocity profile, assuming that the motion is isotropic. This estimated kinetic energy is $`3\times 10^{51}`$ ergs. We can also determine the H I kinetic energy using the total H I mass and the expansion velocity of the superbubble shell; the resulting kinetic energy is $`3.7\times 10^{51}`$ ergs, compatible with the H I kinetic energy estimated from the second moment of the velocity profile. Note that these H I mass and kinetic energy estimates should be considered as upper limits, as some of the H I gas may not be physically associated with the superbubble.

## 4 Stellar Energy Injected into DEM L 192

Stellar mechanical energy is injected into the ambient ISM via fast stellar winds and supernova explosions. Below we determine the contributions of these two types of stellar energy feedback.

### 4.1 Stellar Wind Energy

To determine the stellar wind energy injected into DEM L 192, we consider the 50 most luminous stars in the OB associations LH 51 and LH 54 given by Oey & Smedley (1998). The spectroscopically determined spectral types of these stars are matched with the theoretical models of Schaerer & de Koter (1997) to assess the stellar mass $`M`$, mass loss rate $`\dot{M}`$, and stellar wind terminal velocity $`v_{\mathrm{}}`$. Since $`Z_{\mathrm{LMC}}1/3Z_{}`$, we apply the necessary corrections to the listed theoretical values of $`\dot{M}`$ and $`v_{\mathrm{}}`$ for low metallicity stars as described in Schaerer & de Koter (1997). From these data, we estimate the stellar wind luminosity by $`L_\mathrm{w}=(1/2)\dot{M}v_{\mathrm{}}^2`$. The metallicity-corrected parameters and mechanical luminosities of the 50 stars are listed in Table 1 in the same order as given in Table 1 of Oey & Smedley (1998). We calculate the total mass loss rate to be $`4.4\times 10^5`$ M yr-1 and the total stellar wind luminosity to be $`5.5\times 10^{37}`$ ergs s-1. Thus, for an age of $`3`$ Myr (Oey & Smedley, 1998), the total stellar mass loss is $`100`$ M, and the total wind energy input is about (5 $`\pm `$ 1) $`\times `$ 1051 ergs.

### 4.2 Supernova Energy

To estimate the number of stars that have already exploded as supernovae, we use the combined mass function of LH 51 and LH 54. $`UBV`$ photometry of 456 stars in LH 51 and LH 54 is available from Oey (1996a). The mass of each main sequence star is determined by comparing its location on a color-magnitude diagram (CMD) to the evolutionary tracks of stars with different masses. As the reddening is not known for individual stars, we use the extinction-free magnitude $`W=VDM(BV)A_V/E(BV)`$ and reddening-free color $`Q=(UB)(BV)E(UB)/E(BV)`$, where $`DM`$ = 18.5 is the distance modulus of the LMC, $`E(UB)/E(BV)=0.72`$, and $`A_V/E(BV)=3.1`$. Fig. 4 shows the reddening-free CMD of LH 51 and LH 54, where the evolutionary tracks of main sequence stars with $`Z=0.4Z_{}`$ provided by Lejeune & Schaerer (2001) are plotted for initial masses of 5, 7, 10, 12, 15, 20, 25, 40, and 60 M. Stars in each mass bin are counted. We assume a Salpeter initial mass function for LH 51 and LH 54: $`f(M)=KM^{2.35}`$, where $`K`$ is a constant that can be determined from the star counts. The number of stars with masses between Mmin and Mmax is $`N=_{\mathrm{M}_{\mathrm{min}}}^{\mathrm{M}_{\mathrm{max}}}f(M)dM`$. For a Mmin of 5–7 M and a Mmax of 20–25 M, we find K = 3040 $`\pm `$ 460. This constant is insensitive to the choice of Mmax but rather sensitive to the choice of Mmin; the error bar reflects the uncertainty caused by different Mmin. The number of stars with masses greater than 25 M is $`_{25\mathrm{M}_{}}^{100\mathrm{M}_{}}f(M)dM=25\pm 4`$. From Table 1, we find 12 stars with masses greater than 25 M, so the number of stars that have exploded in supernovae is $`(25\pm 4)12=13\pm 4`$. We have chosen to count main-sequence stars in the 5–25 M range because the evolutionary tracks of such stars are well constrained by the empirically derived relationship between stellar masses and spectral types and are not subject to the uncertainties for the most massive stars (e.g., Aloisi et al., 1999). We have compared the evolutionary tracks of stars in the 5–25 M range between the Geneva models (Schaerer, Meynet, Maeder, & Schaller, 1993) and the Padua models (Fagotto et al., 1994), and find that they are indeed similar. Assuming that each supernova releases about 1051 ergs of explosion energy, the total supernova energy input is approximately (1.3 $`\pm `$ 0.4) $`\times `$ 1052 ergs.

## 5 Energy Budget of the Superbubble DEM L 192

The energies stored in the superbubble DEM L 192 and the stellar energy input are summarized in Table 2. The values and error bars of the thermal and kinetic energies take into account the most probable range of filling factors. The stellar energy input includes contributions from fast stellar winds and supernovae. The superbubble energy includes the thermal energy of its hot interior and the kinetic energy of the superbubble shell.333The thermal energy of the ionized superbubble shell is not included in the energy budget discussion because the superbubble shell is photoionized and the thermal energy is provided by the UV radiation of the OB associations. It is immediately clear that the expected amount of stellar energy injected into the ISM is $`3`$ times as large as that observed in the superbubble.

Previously, observed dynamics of bubbles blown by single stars or superbubbles have been compared with pressure-driven bubble models by Weaver et al. (1977) or variations of such models, and it has been commonly concluded that the observed stellar wind strength is too high for the observed bubble dynamics or that the bubble/superbubble shell expands too slowly for the given stellar wind strengths (e.g., García-Segura, Langer, & Mac Low, 1996; Oey, 1996b; Nazé et al., 2001). These results are equivalent to what we have found in DEM L 192, i.e., the stellar mechanical energy injected into the ISM exceeds the energies storied in the superbubble. More specifically, according to the Weaver et al. (1977) model, the 35 km s-1 expansion velocity and 60 pc radius of DEM L 192 indicate a dynamic age of 1 Myr, which is shorter than the age of the OB associations LH 51 and LH 54. The mass of the ionized superbubble shell consists of the material swept up from the superbubble interior, thus the volume and mass of the H II gas given in §3.2 imply an ambient density of $`3`$ H-atom cm-3. The wind mechanical luminosity required by Weaver et al. (1977) for DEM L 192 is $`1.5\times 10^{38}`$ ergs s-1, and the total energy required during the 1 Myr dynamic time is $`5\times 10^{51}`$ ergs, which is much smaller than the total stellar energy expected from LH 51 and LH 54.

Evidently there is an energy crisis: the observed interstellar energies cannot account for all the stellar energy injected into the ISM. Several mechanisms may drain energies from a superbubble. The most obvious mechanism is a superbubble blowout or breakout, in which the interior hot gas spews out, directly reducing the thermal energy and pressure of the superbubble interior. However, no blowouts are obvious in DEM L 192. As shown in Fig. 1, the extended diffuse X-ray emission from DEM L 192 is completely bounded by outer, curved H$`\alpha `$ filaments, in sharp contrast to the blowout of the superbubble DEM L 152, also known as N44, in which the diffuse X-ray emission extends beyond open, streamer-like H$`\alpha `$ filaments (Chu et al., 1993).

Another mechanism that may drain the internal energy of a superbubble is the evaporation of dense cooler gas into the hot interior. In addition to the dense ionized superbubble shell, DEM L 192 contains dense cloudlets with ionized surface layers, as revealed in high-resolution HST WFPC2 H$`\alpha `$ images (Chen et al., 2000). The main energy loss would occur at the interface layer between the $`10^6`$ K hot gas and the $`10^4`$ K cool gas because the cooling rate is the highest for temperatures of a few $`\times 10^5`$ K. We may estimate the energy loss due to cooling from interface layers. Assuming that the interface layers are in pressure equilibrium with the hot interior gas, their electron number density would be 0.3–0.4 cm-3. The total volume of the interface layers is nevertheless unknown because no interface layer has ever been spatially resolved observationally and even the existence of heat conduction and mass evaporation may be questionable (Shelton, 2003). The spatial extent of an interface layer can be estimated from the relative locations of diffuse X-ray emission and the optical shell of a bubble or superbubble. The highest-quality X-ray observations of wind-blown bubbles are the Chandra observation of NGC 6888 and the XMM-Newton observations of S308 (Gruendl et al., 2004; Chu et al., 2003). While NGC 6888 shows no appreciable separation between the outer edge of the diffuse X-ray emission and the outer edge of the optical shell, S308 shows a gap of 90–200′′ (or 0.5 to 1.7 pc), which contains the dense ionized gas shell and the interface layer. The interface layer must be less than 0.5 pc in thickness. Adopting this upper limit of the interface layer thickness for DEM L 192, the volume occupied by the interface layer is no more than 2.5% of the total volume. Using a cooling function $`\mathrm{\Lambda }(10^5\mathrm{K})/N_\mathrm{H}^23\times 10^{22}`$ ergs cm3 s-1 appropriate for the LMC metallicity (Dalgarno & McCray, 1972), the current cooling rate would be $`4\times 10^{37}`$ ergs s-1. If the superbubble has been evolving with the radius $`t^\eta `$, where $`t`$ is the age and $`\eta `$ is a constant (e.g., $`\eta `$ = 0.6 for the Weaver et al. model), then the volume has been growing $`t^{3\eta }`$, and the total energy radiated away in 3 Myr by the interface layer would be $`4\times 10^{51}(3\eta +1)^1`$ ergs. If $`\eta `$ = 0.6, the total radiative loss is $`1.5\times 10^{51}`$ ergs. This generous upper limit for radiative loss is non-negligible but not sufficiently significant in the energy budget. Therefore, radiative loss in the interface layer between dense cold gas and the hot superbubble interior cannot be the cause of the energy discrepancy.

Where did the stellar energy go if it is not converted into the thermal and kinetic energies of the ISM? The interaction between a fast moving plasma, being it tenuous stellar winds or dense supernova ejecta, and the ambient ISM is a complex magnetohydrodynamic problem (Spitzer, 1978). It is perhaps an oversimplification to assume that supernova ejecta and fast stellar winds are thermalized in the hot interior of a superbubble. In fact, the XMM-Newton X-ray observations of DEM L 192 reveal a possible mechanism that may play an important role in the energy budget. Our analysis of the X-ray spectrum of the diffuse emission from DEM L 192 indicates that not all emission originates from a thermal plasma; a power-law component is needed to explain the emission in the 1.0–3.0 keV range. This power-law component cannot originate from unresolved point sources, as the smoothed X-ray image in the 1.1–3.5 keV band displayed in Fig. 5 shows that this hard X-ray emission is diffuse. It is possible that this power-law component originates from the interactions of supernova ejecta and fast stellar winds with the hot interior of the superbubble, but the exact origin is unclear. It is unlikely that this power-law component represents a synchrotron X-ray emission similar to those seen in young supernova remnants such as SN 1006 (Koyama et al., 1995) and Cas A (Vink & Laming, 2003) and the superbubble 30 Dor C (Bamba et al., 2004) for two reasons. First, the X-ray photon index of DEM L 192, $`\mathrm{\Gamma }=1.3\pm 0.2`$, is lower than those observed in the other objects, $`\mathrm{\Gamma }2`$. Second, the synchrotron X-ray emissions from young supernova remnants and 30 Dor C show limb-brightened morphologies because the relativistic particles are accelerated in the shocks, but the power-law X-ray emission from DEM L 192 does not show enhancement along the edge (see Fig. 5). Observations of more superbubbles are needed to determine whether a power-law component is commonly needed to model the spectra of their diffuse X-ray emission. Future theoretical models of superbubbles need to solve the plasma interaction problem realistically instead of making simplifying assumptions about shocks and converting mechanical energy into only kinetic and thermal energies, so that true progress in the understanding of bubbles and superbubbles can be made.

YHC acknowledges the support of NASA grant NAG 5-13076. We thank the referee for prompt reading and providing useful suggestions to improve this paper. We also thank Yaël Nazé for her critical reading of the paper.