text:
# Direct observation of irrotational flow and evidence of superfluidity in a rotating Bose-Einstein condensate

## Abstract

We observed the expansion of vortex-free Bose-condensates after their sudden release from a slowly rotating anisotropic trap. Our results show clear experimental evidence of the irrotational flow expected for a superfluid. The expansion from a rotating trap has strong features associated with the superfluid nature of a Bose-condensate, namely that the condensate cannot at any point be cylindrically symmetric with respect to the axis of rotation since such a wavefunction cannot possess angular momentum. Consequently, an initially rotating condensate expands in a distinctively different way to one released from a static trap. We report measurements of this phenomenon in absorption images of the condensate taken along the direction of the rotation axis.

In their recent theoretical paper, Edwards et al.  calculated the time evolution of the expansion of a slowly rotating BEC. At higher rotation rates, greater than about half the radial trap frequency, the superfluid nature of the condensate has been demonstrated strikingly in the formation of quantised vortices . However they point out that the superfluidity of the condensate, together with the conservation of angular moment, also leads to clearly observable consequences at lower rotation rates, under conditions where no vortices are present. We report here the experimental verification of these predictions by time-of-flight measurements on a condensate released from a rotating anisotropic trap . These data give further evidence for the superfluidity of Bose-condensed gases. Other recent experiments on this important aspect of Bose-condensation include the observation of quantized vortices , the absence of dissipation below a critical velocity  and the quenching of the moment of inertia in the scissors mode .

The solutions of the hydrodynamic equations for superfluids (that also describe the scissors mode ) show that a condensate always has a moment inertia less than that of a rigid body of the same mass distribution . For the special case of cylindrical symmetry, the superfluid has zero moment of inertia about the symmetry axis. Thus a condensate with some angular momentum can never have a circular shape when viewed along the rotation axis, since this would imply the unphysical situation of infinite angular velocity. An initially rotating condensate therefore expands in a different way to one released from static trap. We observed this phenomenon in absorption images of the condensate taken along the direction of the rotational axis \- the projection of the cloud in this direction never becomes circular (an aspect ratio of unity). The evolution of the condensate density distribution is calculated using the same hydrodynamic equations as in , and agrees well with our data. We give a brief review of the underlying theoretical aspects before describing the experimental procedure and results.

In our experiment we have an anisotropic harmonic potential with three angular frequencies $`\omega _x<\omega _y<\omega _z`$. The potential rotates about the z-axis with angular frequency $`\mathrm{\Omega }`$. It is convenient to define two new parameters: the trap deformation, $`\lambda =(\omega _x^2\omega _y^2)/(\omega _x^2+\omega _y^2)`$ and the mean of the frequencies in the plane of rotation defined as $`\omega _m^2=(\omega _x^2+\omega _y^2)/2`$. In the hydrodynamic limit, a condensate rotating in this trap displays a quadrupolar flow pattern, as described in fig.1, and has a wavefunction which looks like:

$$\mathrm{\Psi }(𝐫)=\sqrt{\rho (𝐫)}e^{i\frac{m\nu }{\mathrm{}}xy},$$
(1)
where $`\rho (𝐫)`$ is the condensate density given in Eq.3 below. The quadrupole frequency, $`\nu `$, can be found from the solutions of the following cubic equation:

$$\stackrel{~}{\nu }^3+(12\stackrel{~}{\mathrm{\Omega }}^2)\stackrel{~}{\nu }+\lambda \stackrel{~}{\mathrm{\Omega }}=0,$$
(2)
where we introduced the dimensionless quantities $`\stackrel{~}{\nu }=\nu /\omega _m`$ and $`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }/\omega _m`$.

In a rotating trap, the condensate experiences effectively different trapping frequencies due to its quadrupolar motion and the condensate density can be written as:

$$\rho (𝐫)=\frac{\mu }{g}\left[1\frac{m}{2}(\stackrel{~}{\omega }_x^2x^2+\stackrel{~}{\omega }_y^2y^2+\omega _z^2z^2)\right].$$
(3)
The modified frequencies $`\stackrel{~}{\omega }_x`$ and $`\stackrel{~}{\omega }_y`$ are given by:

$`\stackrel{~}{\omega }_x^2`$ $`=`$ $`\omega _x^2+\nu ^22\nu \mathrm{\Omega }`$ (4)
$`\stackrel{~}{\omega }_y^2`$ $`=`$ $`\omega _y^2\nu ^22\nu \mathrm{\Omega }`$ (5)
where $`\nu `$ is found from Eq.2.

Hence, the aspect ratio of the condensate in the trap is

$$\frac{R_x}{R_y}=\frac{\stackrel{~}{\omega }_y}{\stackrel{~}{\omega }_x}=\sqrt{\frac{\mathrm{\Omega }+\nu }{\mathrm{\Omega }\nu }},$$
(6)
where $`R_x,R_y`$ are the sizes of the condensate in the x and y-direction, respectively. The effective chemical potential and condensate sizes in the trap can also be calculated from Eqs.5.

To calculate the expansion of a condensate when released from the anisotropic harmonic potential, we use the following ansatz of spheroidal form for the condensate density $`\rho (𝐫,t)`$ and the velocity field, $`𝐯(𝐫,t)`$, (as in ):

$`\rho (𝐫,t)`$ $`=`$ $`\rho _0(t)\rho _x(t)x^2\rho _y(t)y^2\rho _z(t)z^2\rho _{xy}(t)xy`$ (7)
$`𝐯(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(v_x(t)x^2+v_y(t)y^2+v_z(t)z^2+v_{xy}(t)xy\right)`$ (8)
Inserting this ansatz into the hydrodynamic equations of superfluids yields a set of nine coupled differential equations for the expansion parameters, which we integrate numerically. The initial conditions for the density are found from Eq.3. Note that the cross-term $`\rho _{xy}(0)`$ is assumed to be zero at the instant when the condensate is released from the trap. During the evolution in time of flight, this term grows and can become comparable in size to $`\rho _x`$, so that the elliptical condensate rotates with respect to its initial release angle. All components of the velocity field are initially zero, except for the cross term, which is the quadrupole velocity field of the rotating condensate, given by $`v_{xy}=2\nu `$ (see Eq.1). The behavior of the condensate is thus completely determined by the above initial conditions and the nine coupled differential equations .

The angle of the condensate and its sizes can be found by diagonalising the quadratic form defined in Eq.7. One finds, in time of flight, that at first the condensate expands along its smaller axes until the aspect ratio has reached a critical value close to spherical. The condensate has a quadrupolar flow pattern and an associated angular momentum. Conservation of this angular momentum does not allow the condensate to become spherical, as a spherical quadrupolar flow pattern has no angular momentum. Thus, when approaching this critical value in the aspect ratio, the condensate rotates quickly and continues to expand along its original long axis. The rapid increase in angular velocity close to the critical aspect ratio arises from the decrease of the moment of inertia when the aspect ratio approaches spherical . However, energy conservation limits how fast the condensate can rotate and thus the aspect ratio cannot become smaller than a certain value. This behaviour is in contrast to that of a thermal gas which can become spherical and preserve angular momentum, since there are no constraints of irrotationality. Thus it is essentially the superfluid nature of the Bose-condensate which prevents it from becoming spherical.

To observe this behaviour we evaporatively cooled 87Rb atoms in a TOP trap to a temperature of $`0.5T/T_c`$, where $`T_c`$ is the critical temperature for condensation. This produced Bose-condensates of $`1.5\times 10^4`$ atoms, in a trapping potential with frequencies $`\omega _x/2\pi =\omega _y/2\pi =124`$ Hz and $`\omega _z/2\pi =350`$ Hz. We then made the trap elliptical with $`\omega _y/\omega _x=1.4`$ (corresponding to $`\lambda =0.32`$) by changing the ratio of the two TOP-field components to $`B_x/B_y=4.2`$. The eccentric trap was rotated by modulating the high frequency sinusoidal TOP signal ($`7`$ kHz frequency) at the low trap rotation frequency $`\mathrm{\Omega }`$, as described in . The condensate was spun up in $`500`$ ms during which time the eccentricity was ramped up from zero to its final value. Then the condensate was left in the rotating trap for another $`500`$ ms before it was released at a well defined point in the trap rotation, from which the angles of the expanding cloud were measured. Figure 2 shows typical absorption images of the condensate, taken along the axis of rotation, after different expansion times. For these pictures the trap was rotating at 28 Hz and at the instant of release the long axis of the cloud was horizontal. The angle and aspect ratio of the cloud was obtained from a 2D parabolic fit to the density distribution. Figure 2b shows the minimum aspect ratio of 1.31, reached after $`4`$ ms of expansion. The angle increases steadily to an asymptotic value of $`55`$ degrees in fig.2d.

It was necessary to go to a high trap eccentricity so that a clear deformation of the cloud could be observed. We investigated the evolution of the condensate in time of flight at two different trap rotation frequencies $`\mathrm{\Omega }/2\pi =20`$ and $`28`$ Hz. In both cases the condensate was released from a trap of $`\omega _x/2\pi =60`$ Hz, $`\omega _y/2\pi =1.4\times 60`$ Hz and $`\omega _z/2\pi =206`$ Hz. Fig.3 shows the calculated change of angle of a condensate in time of flight after being released from a trap rotating at $`\mathrm{\Omega }/2\pi =20`$ Hz (solid line) and $`\mathrm{\Omega }/2\pi =28`$ Hz (dotted line), with the experimental data points superimposed. In both cases the angle evolves in a similar manner and reaches $`45`$ degrees after about $`6`$ ms. After $`18`$ ms the angle assumes an asymptotic value between $`55`$ and $`60`$ degrees.

Figure 4 shows the theoretical prediction for the evolution of the condensate’s aspect ratio in time of flight, with experimentally measured values superimposed, for $`\mathrm{\Omega }/2\pi =28`$ Hz (upper curve), for $`\mathrm{\Omega }/2\pi =20`$ Hz (middle curve) and for release from a nonrotating trap (lower curve). The data clearly demonstrate how the aspect ratio for an initially rotating condensate decreases up to a critical point, which is reached after approximately $`4`$ ms. From that point on it does not continue to expand along its minor axis but the aspect ratio increases again because the condensate cannot become circular under these conditions. However, the condensate released from a static trap has no velocity field which prevents it from becoming circular and the aspect ratio decreases steadily from $`\omega _y/\omega _x`$ to a final value of less than one. Every experimental point displayed is the average of several measurements. For each time of flight we had to refocus our imaging system as the atoms move out of focus under gravity when released from the trap. However incorrect focusing would only result in a more circular image and the minimum value for the measurement of the aspect ratio of rotating condensates was never consistent with unity. There is remarkable agreement between the experimental data and the theoretical predictions. However, we observed a small deviation of the experimental data from the predicted values for the higher rotation frequency ($`\mathrm{\Omega }/2\pi =28`$ Hz, shown in the upper curve). We do not know the origin of that deviation.

Our results show that an expanding vortex-free Bose-condensate with some angular momentum refuses to become circular about the axis of rotation, as predicted by Edwards et al . This provides direct evidence that Bose-condensed gases flow irrotationally as a consequence of their superfluidity. These measurements complement the observations of the scissors mode  where the irrotational flow was deduced from measurements of the oscillation frequencies of a trapped condensate - a superfluid has less moment of inertia than a classical object, of the same size, undergoing rigid body motion hence the superfluid oscillates faster. Other measurements of the flow of a Bose-condensed gas would be of great interest, for example flow through a narrow tube , or array of holes, analogous to the superleak experiments with liquid helium.

We thank M. Edwards, S. Stringari and the members of the Oxford theoretical BEC group for fruitful discussions. We are very grateful to J. Arlt for his contribution to the early stages of this experiment. We acknowledge support from the EPSRC, St. John’s College, Oxford (G.H.), Christ Church College, Oxford (E.H.), Linacre College, Oxford, MCFA and the EC (O.M.).


format: text

text:
# CCD OBSERVATIONS OF THE RR LYRAE VARIABLES IN THE GLOBULAR CLUSTER NGC 5897

## 1 INTRODUCTION

NGC 5897 (C1514-208) is a moderately metal poor globular cluster with \[Fe/H\] = $`1.80`$ (Harris 1996). Color-magnitude diagrams for NGC 5897 have been published by Sandage & Katem (1968, hereafter SK), Ferraro et al. (1992), Sarajedini (1992) and Testa et al. (2001) and all of them show that the horizontal branch is predominantly blue. The cluster has nine known variable stars: seven RR Lyrae (Wehlau 1990, hereafter W90), one red variable (Eggen 1972) and one SX Phe (Wehlau et al. 1996). W90 has shown that the RR Lyrae variables have a very unusual period distribution; three have periods between 0d.797 and 0d.856, two at 0d.420 and 0d.454 and two more at 0d.342 and 0d.349. As a result, the mean period of the three RRab variables is 0d.828. The fact that the periods of the RRab variables in NGC 5897 are longer than 0d.60 is not unusual because Clement et al. (2001) have shown that for all clusters more metal poor than $`1.70`$, the mean period of the RRab stars is greater than 0d.60. What is unusual, however, is that $`<P_{ab}>`$ is greater than 0d.80; the next longest $`<P_{ab}>`$ for a metal poor cluster is 0d.708 for NGC 4833 which has \[Fe/H\]$`=1.79.`$111Generally, the RR Lyrae variables in metal poor clusters have longer periods than the ones in metal rich clusters. However, Pritzl et al. (2000) have recently shown that the two metal rich clusters NGC 6388 and NGC 6441 have $`<P_{ab}>`$ of 0d.71 and 0d.76, respectively, for their RR Lyrae variables. The variables in these two clusters seem to be brighter than solar neighbourhood field RR Lyrae stars of comparable \[Fe/H\], but the reason for this is not yet fully understood.

Another curious feature of NGC 5897 is the fact that there appear to be nonvariable horizontal branch stars located in its RR Lyrae gap. Both SK and W90 have shown that the stars SK 116 and SK 174 have colors and magnitudes comparable to the RRc stars and that SK 120 lies in the instability strip well to the red of the RRc stars and to the blue of the RRab stars. Smith (1985) made a spectroscopic study of SK 120 and concluded that the star is probably a cluster member. We therefore decided it would be worthwhile to observe NGC 5897 with a CCD detector. It may turn out that some of these ‘apparent’ nonvariable stars are in fact varying, but with an amplitude lower than can be readily detected in photographic studies. All of the published photometric studies of the RR Lyrae variables in NGC 5897, with the exception of a set of CCD observations of V2 by Wehlau et al. (1996), are based on photographic data.

## 2 THE OBSERVATIONS

Our investigation is based on 526 CCD frames obtained with the University of Toronto’s 61 cm ($`f/15`$) Helen Sawyer Hogg (HSH) telescope at the Las Campanas Observatory of the Carnegie Institution of Washington. The observations were made on six consecutive nights (May 1 to 6, 1997) with the Kodak KAF-4200 CCD system that Horch et al. (1997) and Slawson et al. (1999) used with the same telescope in February 1997. The chip has 9 $`\mu `$m square pixels arranged in a 2033$`\times `$2048 format, but before readout, the data were binned $`2\times 2`$ so that the effective pixel size for our data was 18 $`\mu `$m square. The field of view (shown in Figure 1) is approximately $`6.9\times 6.9`$ arcminutes and includes all of the RR Lyrae variables studied by W90. The frames were obtained through a V filter and all of the exposure times were 2 minutes. These short exposures were necessary because of severe tracking problems with the telescope.

The frames were cleaned using IRAF. Then, for the photometric reductions, we combined the frames in groups of three, so that we could obtain photometry with a precision of $`0.02`$ mag (or better) at V$`16.5`$, the level of the horizontal branch. In Table 1, we list a sample of the data obtained on the night of May 1 to illustrate how the frames were combined. The photometric reductions were carried out using the IRAF versions of DAOPhot-II and ALLSTAR-II and the stand-alone version of Stetson’s (1987, 1993, 1994) ALLFRAME code. Approximately 1000 stars were identified on each frame.

In order to convert our instrumental magnitudes to standard V magnitudes, we used the photometry of SK. Since none of their photoelectric standards were in the field shown in Figure 1, we obtained 7 frames on the night of May 4 with the telescope offset to a position $`5.7`$ arcminutes further north. With this offset, we were able to include the SK photoelectric standard stars A, C, E, K, M, N, O and P in the field of view. We determined the instrumental magnitudes for the stars that were in the area overlapping the two fields on the same magnitude scale as the ‘offset’ frames and then calculated magnitudes on the same scale for all of the program stars. Finally, the transformation from instrumental to standard V magnitudes was derived for SK’s photoelectic standards and then applied to all of the other stars.

## 3 THE VARIABLES

### 3.1 Search for Variables and Period Determination

The known RR Lyrae variables were identified and found to have a mean $`V`$ magnitude of $`16.25`$. To search for additional RR Lyrae variables, we examined the observations for all of the stars with mean $`V`$ magnitude between $`15.75`$ and $`16.75`$ (approximately 150 stars). We plotted light curves to check for periodic variations. As a result, four new RR Lyrae variables were identified. We have numbered them V10 to V13 and they are labelled in Figure 1. This increases the total number of known RR Lyrae variables in NGC 5897 from seven to eleven. In order to determine the periods for the new variables, we used Stellingwerf’s (1978) phase dispersion minimization technique to estimate an approximate value and then made small adjustments so that the observations for the different nights would be well aligned. Our adopted periods are listed in Table 2. For the known variables, we adopted the periods listed by W90. Also included in the table are the SK number, the x,y position on the frame in pixels, the x,y position relative to the cluster center in arcseconds on the system of Sawyer Hogg (1973), the mean $`V`$ magnitude, the $`V`$ amplitude along with the standard deviation of the fit to equation (1) and the Fourier phase differences $`\varphi _{21}`$, $`\varphi _{31}`$ along with their standard errors. The mean magnitudes and amplitudes were derived from a fit to a Fourier series of the form:

$`mag=A_0+{\displaystyle \underset{j=1,n}{}}A_j\mathrm{cos}(j\omega t+\varphi _j)`$ (1)
where $`\omega `$ is ($`2\pi `$/period). The Fourier phase differences ($`\varphi _{ji}`$) are ($`\varphi _jj\varphi _i`$) from equation (1). For each star, the epoch was taken as HJD $`2,450,500`$ so that $`t`$ in the equation refers to (HJD$`2,450,500`$) and HJD represents the mean heliocentric Julian date of the observation (in this case, the combination of three frames). Light curves, arranged in order of increasing period, are shown in Figure 2. A sample of the observations for V1 is shown in Table 3. All of the observations for V1-3, V6-8 and V10-13 are presented in the electronic version of the table. V4 was not included in our study because its image was blended with that of another star and as a result we were unable to obtain photometry with the desired precision. In fact, many of the frames could not be used to study the other variables either, because of the tracking problems.

It turns out that two of the new variables, V11 and V13 are SK 120 and SK 116 respectively, apparent nonvariable stars that lie in or near the RR Lyrae gap according to SK and W90. Since the $`V`$ amplitudes for these two stars are less than $`0.20`$ mag, they were not detected in photographic studies of the cluster. We also examined the magnitudes for SK 174, a star that W90 found to be near the blue edge of the instability strip, but no systematic variation was detected. Its mean $`V`$ magnitude is $`16.392`$ with a standard deviation of $`0.029`$. In comparison, the standard deviation of the previously known variables ranged from $`0.053`$ for V8 to $`0.215`$ for V1.

The other two new variables, V10 and V12, have larger amplitudes, but their periods are close to simple fractions of a day. This probably accounts for the fact that they were not discovered previously. In fact, even in the present study, the phase coverage for V10 (P$``$0d.5) is incomplete. We considered the possibility that V10 has a period twice as long, but concluded this to be unlikely. The shape of its light curve is similar to that of three $`\omega `$ Centauri variables with similar periods, V47, V68 and V123 (OGLE ID 185, 73 and 169), observed by Kaluzny et al. (1997, hereafter K97). All of these stars have an inflection on the rising branch of their light curves, just before maximum light. Furthermore, the $`\omega `$ Cen stars for which K97 found periods of about 1 day are all more than $`0.5`$ mag brighter than the RR Lyrae variables. Since V10 is only $`0.2`$ magnitudes brighter than the mean $`<V>`$ for the other variables, we assume that our adopted period is the correct one.

### 3.2 Period-Luminosity and Period-Amplitude Relations

The period-luminosity and period-amplitude relations are plotted in Figure 3. In the upper panel of the diagram, the stars separate into three separate regions which we interprete to be due to different modes of pulsation: fundamental, first and second-overtone. W90 has already noted that the variables with periods $``$0d.34 may be pulsating in the second-overtone mode. Another interesting feature of Figure 3 is its similarity to the P-L and P-A plots that Clement & Rowe (2000, hereafter CR) made for the $`\omega `$ Centauri variables observed by K97. To illustrate this, we plot in Figure 4 the Fourier phase differences $`\varphi _{21}`$, $`\varphi _{31}`$ and the $`V`$ amplitude against $`\mathrm{log}P`$ for the NGC 5897 variables and include the $`\omega `$ Cen variables that CR considered to belong to the Oosterhoff (1939, 1944) type II class (i.e. $`<V>14.65`$ for the fundamental mode and $`<V>14.60`$ for the overtone modes). Also included is star 96 ($`\mathrm{log}P=0.3`$, $`A_V=0.17`$) of K97 which was not included in CR’s study.

In Figure 4, it can be readily seen that the RRab stars V1 and V6 have periods, Fourier phase differences and amplitudes comparable to some of the $`\omega `$ Cen Oosterhoff type II fundamental mode pulsators. What is striking though is that V1 and V6 have longer periods and hence lower amplitudes than most of the $`\omega `$ Cen stars. A similar trend can be seen among the first-overtone variables. Both V2 and V11, with periods $``$0d.46 ($`\mathrm{log}P0.34`$), have much lower amplitudes than most of the $`\omega `$ Cen first-overtone pulsators. The case of V2 is particularly interesting because W90 found an increase in period and a decrease in amplitude for this star between the 1950s and the 1980s. The $`V`$ amplitude we have derived for V2 ($`0.24`$ mag) is comparable to her lower, more recent value ($`0.22`$ mag), but the $`V`$ amplitude that she derived for the 1956-1966 observations was $`0.39`$ mag, similar to the amplitudes of the $`\omega `$ Cen first-overtone variables. V2 seems to have changed significantly during the last 50 years. We analysed our observations to search for evidence that V2 is in the process of mode switching. To do this, we measured the residuals to the light curve and searched for periods in the range between 0d.25 and 0d.65. No oscillations with an amplitude greater than $`0.03`$ mag were detected for any period in the interval. If the star were in the process of changing from first-overtone to fundamental mode pulsation, we would have expected to detect oscillations with a period $``$0d.60. We therefore conclude that if V2 is in the process of switching modes, the fundamental mode oscillations are still very weak.

The four stars with P$``$0d.34 (V7, V8, V12 and V13) have periods and amplitudes comparable to some of the $`\omega `$ Cen second-overtone variables, but, like the RRab stars, they seem to congregate at the long period end of the sequence in the P-A plot. CR found that several of the second-overtone variables in $`\omega `$ Cen exhibited non-radial pulsations. We performed a period search on the residuals to the light curves of the above-mentioned four stars in an attempt to search for non-radial pulsations. No oscillations with amplitudes greater than $`0.03`$ mag were detected in any of them. However, the amplitude of the non-radial pulsation detected by CR for star 186 in $`\omega `$ Cen was $`0.10`$ mag, much higher than this.

As noted previously, the light curve for V10 has a structure similar to the light curves of V47, V68 and V123 in $`\omega `$ Cen and we suggest that it might be an anomalous Cepheid. The properties of anomalous Cepheids (ACs) have been discussed by Nemec et al. (1994, hereafter N94). They are located in the Cepheid instability strip and are more luminous than horizontal branch stars, but they tend to have shorter periods than globular cluster Cepheids because they are more massive. They are probably coalesced binary stars and may be related to blue stragglers. Most of the known ACs occur in nearby dwarf galaxies, but N94 pointed out that V68 in $`\omega `$ Cen might be an AC pulsating in the fundamental mode.222N94 also considered that V84 in $`\omega `$ Cen might be an AC candidate, but in the meantime, van Leeuwen et al. (2000) have shown that V84 is not a member of $`\omega `$ Cen. We assume that if V68 is an AC, then V47 and V123 in $`\omega `$ Cen should also be classified as ACs and if clusters with blue stragglers are more likely to contain ACs, then it is plausible that there are ACs in NGC 5897 as well. Testa et al. (2001) pointed out that the cluster has a sparse, but clearly visible population of blue stragglers. If V10 is an AC, then V3 could be one as well. N94 derived a slope for the period luminosity relation for AC variables: $`\mathrm{\Delta }V/\mathrm{\Delta }\mathrm{log}P=3.13\pm 0.28`$ and if both V3 and V10 are ACs, V3 should be $`0.23`$ mag fainter than V10. The actual difference is $`0.27`$ mag, in reasonable agreement with the prediction of N94. Another possibility raised by N94 is that stars like this might be first-overtone population II Cepheids. Thus, there is still some doubt about the true nature of V3 and V10.

Our study reinforces the conclusion that the period distribution of the RR Lyrae variables in NGC 5897 is unusual. W90 plotted a period-frequency distribution that demonstrated that the cluster has no variables with periods between 0d.45 and 0d.75 where the majority of RRab variables are found in other clusters. Although we have discovered four new variables in this investigation, and two of them have periods between 0d.45 and 0d.50, neither is an RRab star. Furthermore, there are still no known variables with periods between 0d.50 and 0d.75 and most cluster RRab stars have periods in this range.

We thank Adam Muzzin for his assistance in the preparation of the paper. Support for this work from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged.


text:
# Tidal Torques and Galactic Warps

## 1 Introduction

Many edge-on disk galaxies show integral-sign warps, where the majority of the disk is planar but where the outer region of the disk lies above the plane on one side of the galaxy and below the plane on the other (e.g. Binney 1992). Most extended HI disks appear warped (e.g. Briggs 1990) and half of all disk galaxies have optical warps (Reshetnikov & Combes 1998).

Various methods have been proposed for creating and maintaining warps, such as normal bending modes (e.g. Sparke & Casertano 1988) and disks askew in flattened dark matter halos (Toomre 1983; Dekel & Shlosman 1983). Motivated by the idea that infalling material will alter the direction of the angular momentum of a galaxy (Quinn & Binney 1992), Ostriker & Binney (1989) studied the reaction of massive rings to a slewing disk potential. They found that warps occurred in regions of low surface density. Also motivated by the cosmic infall of angular momentum, Debattista & Sellwood (1999) found that when the angular momenta of a halo and disk are misaligned, dynamical friction between them can produce a warp.

In a cosmological setting, a galactic disk is expected to continuously experience tidal torques from a variety of sources, the three more important being the distribution of mass in the local environment of a galaxy, substructure in the dark matter halo such as dwarf satellites and high velocity clouds (HVCs), and a misalignment between the disk angular momentum and the figure axes of the dark matter halo. In this paper, we use cosmological N-body simulations to deduce what gravitational tidal torques a typical galaxy experiences from these three sources, and study whether these torques provide a possible origin for warped disks.

## 2 Torques in a Cosmological Context

We determine the strength and radial dependence of the torque by analyzing a sample of 19 halos formed in an N-body simulation of structure formation in the standard cold dark matter scenario ($`\mathrm{\Omega }_0=1`$; for details, see Navarro & Steinmetz 1997). Each halo consists of at least 1000 particles and has been selected not to be significantly disturbed by recent merger events. We assume a disk to form in the plane perpendicular to the angular momentum of the dark matter halo.

The magnitude and radial dependence of the torque from all sources is determined by calculating the gravitational force component that is exerted by all particles in the simulation and that is acting perpendicular to the disk. The torques are expressed in terms of the tilting timescale $`t_{\mathrm{tilt}}(r)`$ defined as the time required for a solid massless ring of radius $`r`$ to react to the torque by tilting by one radian.

We separate the mass in the simulation into that inside the halo (defined by the virial radius $`r_{200}`$) and that outside the halo and compare the torques due to each. We find the torque due to the mass outside the halo to be negligible, indicating that galaxies in the local universe and cosmic tidal shear do not contribute to the torque a (field) disk galaxy experiences.

The torques inside three sample halos are shown in Figure 1. The solid line is a power-law fit for each halo. The fits reproduce the broad behavior of the torquing forces, but there is significant deviation at specific radii. This power-law type radial dependence is characteristic of misalignments typical for dark matter halos that form by hierarchical clustering (Cole & Lacey 1996; Warren et al. 1992), with mild misalignments of typically 20 for a moderately flattened ($`b/a0.8`$) dark matter halo. The deviations in the radial behavior from the power-law type dependence are due to substructure within the halo, such as satellite dwarfs, and, to some extent, due to numerical noise.

The inner regions of the halo exhibit stronger torques (or equivalently shorter tilting timescales) than the outer parts, i.e., the inner disk should tilt faster. This is because the torquing force, which scales roughly as the density, increases more rapidly toward the centre than the disk’s ability to resist the torque due to its angular momentum. Also note that the timescales are less than a Hubble time (indicated by the dotted lines) over much of the disk, so we should see the effects of these torques in a large number of observed galaxies, consistent with observational findings (Briggs 1990, Reshetnikov & Combes 1998).

## 3 Reaction of the Disk

We study the reaction of a disk to these torques by performing numerical N-body simulations of a massive Milky Way type galactic disk (maximum circular velocity $`v_{\mathrm{max}}=233`$km/s, scale length $`r_d=3.5`$kpc, vertical scale height $`h_z=325`$pc) subject to the torques found in section 2. Equilibrium disk models of disk mass $`1\times 10^{10}M_{}`$, $`3\times 10^{10}M_{}`$, and $`5.6\times 10^{10}M_{}`$ in a static spherically-symmetric NFW halo potential (Navarro, Frenk & White 1997) with concentration parameter $`c=15`$ and a virial velocity of $`v_{200}=175`$km/s are constructed using the method of Hernquist (1993). Each disk contains 16384 particles. The models were evolved using GRAPESPH (Steinmetz 1996) for 2 Gyr, which took 5000–7000 timesteps depending on the model. A plummer softening of 0.3 kpc has been used.

Figure 2a shows the simulation of a $`3\times 10^{10}M_{}`$ disk subjected to a torque for 1 Gyr. The disk was aligned with the $`xy`$-plane at $`t=0`$. The inner part of the disk is flat and clearly tilted toward the positive $`x`$-axis. Beyond 11 kpc, the disk warps back toward the original plane. It distinctly resembles observed warped galaxies. The particles of the simulated galaxy were binned into spherical shells 2 kpc wide, and the minor axis of each bin was found from the moment of inertia tensor. Figure 2b plots the tilt angle of each ring from the initial plane of the disk. The flat region shows up clearly as the inner rings which are all tilted a uniform 17 from the initial plane, while beyond 11 kpc (the warp radius) the disk warps back toward the initial plane. A massless disk, in contrast, does not exhibit an inner flat region and is warped at all radii, as shown by the solid line in Figure 2b. The self-gravity of the massive disk maintains its flatness in its inner regions.

The warp radius $`r_w`$ at which the warp starts, 11 kpc in the disk shown in Figure 2, is well-defined and can be followed over time. It coincides with the radius at which the difference between the tilt of the massive disk and the tilt of an equivalent massless disk is maximized, i.e. the data point that lies highest above the solid line in Figure 2b. As the disk evolves under the influence of the torque, this warp radius moves out through the disk at a rate that depends on the mass of the disk, as seen in Figure 3a. The three sets of symbols represent simulations with disks of different mass. Warps move faster through higher mass disks than through lower mass disks.

The effect of the disk mass on the rate of warp growth implies that the self-gravity of the disk is important to the formation of the warp. In a massive disk, particles at different radii are coupled to each other gravitationally and act to keep the disk locally flat. This suggests that the crucial factor that determines where the warp develops is the local surface density of the disk. In fact, Ostriker & Binney (1989) examined the effect of a slewing disk potential on a set of self-gravitating rings and found that regions of high surface density react like a solid body, but that warps can occur where the surface density is lower. Hofner & Sparke (1994) noted that in most cases the group speed of bending waves in a disk of surface density $`\mathrm{\Sigma }(r)`$ and angular rotation velocity $`\mathrm{\Omega }(r)`$ is $`c_g=\frac{\pi G\mathrm{\Sigma }(r)}{\mathrm{\Omega }(r)}`$, and therefore the time for a warp to settle at a given radius is inversely proportional to the surface density at that radius.

To test this hypothesis, we translate the warp radii $`r_w`$ of Figure 3a into local surface densities using $`\mathrm{\Sigma }(r_w)=\frac{M_d}{4\pi r_e^2}e^{r_w/r_e}`$ for a disk of mass $`M_d`$ and exponential scale length $`r_e`$, and plot the surface density at the warp radius as a function of time for the three disks of different mass in Figure 3b. The surface density at the warp radius falls as the warp moves out through the disk. Note that the local surface densities at the warp *at a given time* are quite similar for all models. It appears that the local surface density is the important parameter for determining the warp radius at a given time.

## 4 Summary

Our conclusions are the following:

* Cosmological N-body simulations show that galactic disks are usually misaligned with the mass distribution of the dark matter halo in which they are embedded. Tidal torques are dominated by this misalignment and by asymmetries in the halo mass distribution, such as satellite dwarfs. The resulting radial dependence of the torque follows a power law of index -2.5. The gravitational torque is typically strong enough to tilt the disk by 1 radian in 1–10 Gyr. Torques from mass outside the halo appear unimportant.
* The tilting of the inner portions of the disk owing to the torque proceeds faster than that of the outer regions, resulting in a warped disk. The self-gravity keeps the inner regions of a massive disk locally flat.
* The radius inside which the disk is flat grows with time. More massive disks have faster-growing warps, mostly because only at large radius is the surface density sufficiently low that its self-gravity cannot maintain the flat disk.

This work was supported by NSF grants AST98-07151 and PHY99-0749, and NSERC grant PGSB-233028-2000.

References

Binney, J. 1992, ARA&A 30, 51

Briggs, F. 1990, ApJ 352, 15

Cole, S., Lacey, C. 1996, MNRAS 281, 716

Debattista, V.P., Sellwood, J.A. 1999, ApJL 513, 107

Dekel, A., Shlosman, I. 1983, in IAU Symposium 100, Internal Kinematics and Dynamics of Galaxies, ed. Athanassoula, E. Dordrecht: Kluwer, 187

Hernquist, L. 1993, ApJS 86, 389

Hofner, P., Sparke, L.S. 1994, ApJ 428, 466

Navarro, J.F., Frenk, C.S., White, S.D.M. 1997, ApJ 490, 493

Navarro, J.F., Steinmetz, M. 1997, ApJ 478, 13

Ostriker, E.C., Binney, J.J. 1989, MNRAS 237, 785

Quinn, T., Binney, J. 1992, MNRAS 255, 729

Reshetnikov, V., Combes, F. 1998, A&A 337, 9

Sparke, L.S., Casertano, S. 1988, MNRAS 234, 873

Steinmetz, M. 1996, MNRAS 278, 1005

Toomre, A. 1983, in IAU Symposium 100, Internal Kinematics and Dynamics of Galaxies, ed. Athanassoula, E. Dordrecht: Kluwer, 177

Warren, M.S., Quinn, P.J., Salmon, J.k., Zurek, W.H. 1992, ApJ 399, 405


text:
# Comment on “Critical current density from magnetization hysteresis data using the critical-state model”

## Abstract

A recent paper \[Physical Review B 64 014508 (2001)\] claims to present an exact method to extract the critical current density JC from the M-H hysteresis curve. We show that this claim appears unjustified.

Keywords: Critical state model, critical current density, M-H hysteresis.

The hysteresis in the isothermal M-H curves of superconductors was first related to the critical current density JC by Bean’s critical state model . Fietz and Webb  showed that the errors involved in extracting the field-dependent JC(B) from these M-H curves are of order (d2JC/dB2). The applicability of this method of extracting JC(B) has been studied in great detail  after the discovery of high-TC superconductors. The recent claim  of an ‘exact method’ of obtaining the same information from the same data, consequently, needs to be examined.

The ‘exact method’ of reference relates, in its equations (7) and (8), the derivatives of the forward and reverse envelope M-H curves to JC(Ba), and to the difference between the applied field Ba and the field BC+ (or BC- ) at the centre of the sample during the increasing (or decreasing) field cycle. It then presents a method (see figures 2 and 3 of reference) to relate the central fields BC+ (or BC-) to Ba via the measured M-H curves. All these results are obtained for a sample in the shape of an infinite slab parallel to the applied field.

Analytical solutions of Bean’s critical state model, for various JC(B), are known for infinite cylinders (of various cross-sections) in parallel field geometry, when the sample demagnetization factor N=0 . It is known that while field profiles B(x) depend on JC(B) and not on the shape of the cylinder’s cross-section, the value of the magnetization M is dictated by the cylinder’s cross-section . While equations (5) and (6) of reference are valid for all N=0 cylinders, the derivation leading to equations (7) and (8) breaks down. This is simply because the volume element at x depends on the shape of the cylinder’s cross-section, resulting in different x-dependent weight-functions within the integral sign. Equations(11) of reference give explicit analytic expressions for the M-H curves of circular cylinders when JC(B) follows the exponential model JC(B) = JC(0)exp(-B/H0). It is also known  for the exponential model that BC+ = H0 log\[exp(Ba/H0) - (H/H0)\] and that BC- = H0 log\[exp(Ba/H0) + (H/H0)\]. It is straightforward to confirm that equations (7) and (8) of reference are not satisfied.

The argument leading to figure 3(b) of reference is based on the relation (depicted in figure 2 thereof) that the total flux contained in the sample is the same in the field-increasing case when the applied field is Ba, as in the field-decreasing case when the applied field is BC+ (corresponding to that Ba). It follows from the results in reference that this is not true for circular cylinders for JC(B) following the exponential or the Kim-Anderson models. It follows from reference that this is not true for a circular cylinder even when JC(B) is independent of B. It is straightforward to see that this result also breaks down for cylinders with other cross-sections .

We assert that the method of reference is applicable only for an infinite slab in parallel field, and is not valid for any other sample shape. The conventional method  for extracting JC(B) from isothermal M-H hysteresis curves is valid for various shapes, and involves errors whose upper bounds are known.

Finally, reference claims to show that the shape of the hysteresis bubble identifying the peak-effect in JC(B) is asymmetric about the peak field BP. This is shown using the input that JC(B) decays as the magnitude of (B-BP) rises. While this result has been shown earlier explicitly with reference to the peak effect (see figure 2b of reference), the calculation is identical to the special case of BP=0. This asymmetry in the ‘central peak position’ has been studied extensively in recent years .


text:
# Hinged Kite Mirror Dissection

## 1 Introduction

A dissection of one polygon to another is a partition of the first polygon into smaller polygonal pieces that can be rearranged to form the second polygon. Dissections are possible when (and only when) the two polygons have the same area, indeed dissection was used by Hilbert as the basis for an axiomatization of area . Dissection puzzles are also popular in recreational geometry, where the main aim is to find dissections of interesting shapes such as regular polygons that use as few pieces as possible .

A classic example is the four-piece dissection from an equilateral triangle to a square, often ascribed to Dudeney (but see  for speculation on its origin). The four pieces in the dissection can be connected by hinges, points of attachment at which the two attached pieces are free to rotate (Figure 1), while still allowing all four pieces to fold up into both the square and the triangle. This example has sparked much interest in similar hinged dissections  but few general results are known, and it remains open whether each equal-area pair of polygons has a hinged dissection. In one of the few theoretical papers in this area, Demaine et al.  showed that chains of isosceles right triangles form hinged dissections between any pair of $`n`$-ominos (Figure 2), and more generally that hinged dissections are possibly between many pairs of polyforms, shapes formed by face-to-face gluing of many copies of the same basic form.

In this paper we demonstrate the existence of another class of hinged dissections: we prove that any asymmetric polygon (including polygons with holes) has a hinged dissection to its mirror image. Our method is to find a dissection that can be unfolded on its hinges into a symmetric form: a chain of kites (quadrilaterals with reflection symmetry across a diagonal) connected end-to-end, so all their axes of symmetry lie on a common line. Clearly, one can then perform the mirror image of the unfolding process to fold the chain back up into the mirror image of the original polygon; Figure 3 shows this unfolding and refolding process applied to a hinged kite dissection of an asymmetric concave hexagon (however this dissection is not one that would be found by the algorithms we describe). We prove that any polygon with $`n`$ sides has a hinged kite dissection with $`O(n)`$ pieces, which can be computed by an algorithm running in $`O(nlogn)`$ time.

Although the set of dissections we find is very limited, we discuss possible implications for more general dissections between any two equal-area polygons. Our technique allows us to reduce the general dissection problem into the question of whether any two equal-area triangles have a hinged dissection in which some copies of a specified pair of vertices of the first triangle map to a specified pair of vertices of the second.

As in the paper of Demaine et al. , we do not consider the question of whether our dissections can be continuously unfolded without any intermediate self-intersections, so in the terminology of Frederickson  all our dissections are wobbly-hinged.

## 2 Dissection Process

We now describe the steps by which we find a hinged dissection of an arbitrary polygon into a chain of kites. These steps are also illustrated in Figure 4. Our technique is based on a circle-packing algorithm of Bern and the author  for partitioning polygons into well-shaped quadrilaterals in the context of finite element mesh generation; this method is based on previous circle-packing nonobtuse triangulation algorithms by Bern et al.  and has also been applied to problems of paper folding and cutting .

We dissect the given polygon by the following sequence of steps. Steps 1-5 are taken from the kite meshing algorithm of Bern and the author, while the remaining steps transform the kite mesh into a hinged dissection. In steps 1-3 we pack the polygon by tangent circles, so that the polygon is partitioned by the circles into regions of two types: interiors of circles, and nonconvex gaps exterior to the circles. Each gap is bounded by three or more sides consisting of segments of polygon boundary and arcs of tangent circles.

1. We begin by placing pairs of congruent circles near each reflex vertex of the polygon, tangent to each other and to the polygon.
2. We place additional circles tangent to each boundary component of the polygon, so that the circles are connected in a cycle by tangencies, with a circle doubly tangent to the polygon near each convex vertex. After this step, each gap involving a segment of polygon boundary (other than the four-sided gaps created in step 1) has exactly three sides. However the gaps in the interior of the polygon may still have many sides.
3. As long as there is a gap with five or more sides, we place a circle to split it into two simpler gaps.
4. The remaining gaps have three or four sides. We draw line segments between each circle center and the circle’s points of tangencies, partitioning the polygon into triangles and quadrilaterals surrounding each gap, with distinguished points (the tangencies) towards the center of each triangle or quadrilateral edge.
5. We now partition each of these triangles or quadrilaterals into kites, according to a case analysis shown in Figure 5:

1. In a three-sided gap interior to the polygon, we place a point at the circumcenter of the triangle formed by the three points of tangency, and connect this center point to each tangency.
    2. In a four-sided gap, the four points of tangency are always cocircular . In most such cases, as in the case of three-sided gaps, we place a point at the circumcenter of these four points, and connect this center point to each tangency.
    3. There may be some four-sided gaps in which the center point is not interior to the convex hull of the four tangencies, so that the previous case would lead to the creation of a concave dart shape instead of a kite. Bern et al.  call this case a bad gap and show that it can always be split into two good gaps by the addition of a single circle tangent to two of the four arcs of the bad gap. These two good gaps can be covered by seven kites (Figure 5, top right).
    4. When two circles form a gap with a straight piece of polygon boundary, we can partition this gap into two kites by a line segment through the circle tangency and perpendicular to the line between the two circle centers. The same type of partition also applies to a gap containing a reflex vertex, because we chose the two circles forming this gap to be congruent.
    5. The final case consists of a gap formed by a convex vertex and a single circle, however this type of gap is already in the form of a kite.
6. We now have a partition of the polygon into kites, however we are not finished because it may not be possible to hinge the kites appropriately. We call the kites of this partition large kites to distinguish them from the ones formed in step 7 below. We next find a tree, with one vertex interior to each large kite, where each tree edge connects points from two adjacent large kites. (I.e., this is a spanning tree of the dual graph of the large kite mesh.)
7. We partition each large kite into four smaller pieces by placing a point at the intersection of its two diagonals, and connecting that point to the midpoints of the large kite edges. This partitions the large kite into four pieces, two of which are similar to the original large kite (shown shaded in Figure 4) and the other two of which are rhombi.
8. We arrange the spanning tree of step 6 so that its vertices lie at the interior points added within each large kite, and its edges lie along the connections from these interior points to the large kite edge midpoints. We add a single segment connecting this spanning tree to the midpoint of an edge on the outer boundary of the polygon.
9. Finally, we trace around the boundary of the tree, and form a linear sequence of small kites and rhombs in the order in which they are visited by this trace. We hinge these kites and rhombs at the vertices on the edge midpoints of the large kites. Each small kite is hinged at the two vertices of its axis of symmetry, and each rhomb is hinged at two opposite vertices.

Bern and the author  use a somewhat more complicated case analysis in step 5, allowing four-sided gaps involving the edges of the polygon, in order to show that the dissection into kites used here can be performed in time $`O(nlogn)`$ and that it need only create $`O(n)`$ pieces.

Another example of a dissection created by this process, of a scalene triangle, is shown in Figure 6 (center). In this case the circle packing consists of a single circle inscribed in the triangle, eventually resulting in a twelve-piece dissection. However, the number of pieces can be improved: as shown in the right of the figure, any scalene triangle has a three-piece hinged mirror dissection into a kite and two isosceles triangles, formed by cutting from the midpoints of the two short sides of the triangle to a third point on the hypotenuse. The third point is the reflection of the hypotenuse midpoint across the perpendicular bisector of the two other midpoints; the line through it and the opposite vertex is perpendicular to the hypotenuse. The three pieces formed by these two cuts are then hinged at the side midpoints.

## 3 Possible Implications

We still seem to be a long way from solving the question of whether hinged dissections exist between any pair of equal area polygons, or even more generally between any set of equal area polygons. However, our kite dissection can be used to reduce this problem to a seemingly more simple form.

Suppose we have two equal area polygons, both dissected into chains of kites hinged end-to-end. The sequences of areas of the kites can be viewed as partitions of the one-dimensional interval $`[\mathrm{0},A]`$, and we can find a common refinement of these two partitions by overlaying them. Geometrically, as shown in Figure 7, this corresponds to introducing a sequence of cuts to the two chains of kites, partitioning them into smaller kites and darts, still hinged end-to-end, so that both chains are composed of polygons that form the same sequence of areas. In other words, the first kite or dart from the first chain has the same area as the first kite or dart from the second chain, and so on. The cuts in one chain correspond to the hinges in the other chain and vice versa. If we could then further hinge-dissect each equal-area pair, we could combine these parts into a hinge dissection of the original two polygons. By further splitting the kites and darts along their axes of symmetry, we reduce the problem to one in which we must dissect a sequence of equal area triangles.

Summarizing, we would be able to hinge-dissect any two equal-area polygons, if only we could hinge-dissect the very simple special case of two equal-area triangles, with the restriction that copies of two vertices from the first triangle are mapped to two vertices of the second triangle so that the dissected triangles can be connected to their neighbors in the chain. Our kite dissection method transforms any single dissection problem of two polygons with a total of $`n`$ sides into a sequence of $`O(n)`$ triangle dissection problems. More generally, we can use the same construction to reduce any $`k`$-way dissection problem to one involving only triangles.

A similar result could be obtained by using Saalfeld’s decomposition of equal-area polygons into combinatorially equivalent equal-area triangulations , however his method lacks complexity bounds and seems to use a large number of pieces.

## 4 Discussion

We have shown that any polygon has a hinged dissection in the form of a chain of kites, that can be unfolded and refolded to form the mirror image of the original polygon. The result also has some possible consequences for the open problem of the existence of hinged dissections between any pair of equal area polygons.

Some questions about our method remain unanswered, for instance whether our dissections or modifications of them can be unfolded in a continuous motion that avoids self-intersections. Also, the number of pieces used by our dissections, although asymptotically optimal, seems large, and Figure 6 shows that it can be reduced by a factor of four in the case of scalene triangles. Is a similar reduction possible more generally?

## Acknowledgements

My thanks go to Greg Frederickson for encouraging me to publish these results, to Erik Demaine for extensive comments on a draft of this paper, and to Cinderella for help with the figures. This work was supported in part by NSF grant CCR-9912338.


text:
# Phase-transitions in a model for the formation of herpes simplex ulcers

## 1 Introduction

One of the most common and intensively studied diseases among humans is the herpes simplex virus (HSV) infection. Apparently, the unique symbiosis that exists in nature between humans and the HSV allows the viral particles to remain inactive (latent infection) in the cranial nerve ganglia after a primary infection, producing frequently recurring localized infections during the host’s lifetime . The reactivation of HSV from latency may occur at any time and it is characterized by active viral replication in the epithelium causing vesicular eruptions in human mucosae and skin. The rupture of these vesicles and the consequent cell necrosis leave the characteristic herpetic lesion or ulcer.

Basically, there are two distinct types of herpes simplex virus, namely, HSV Type I and HSV Type II. The former generally involves infection above the waist (ocular and facial) while the latter infects tissues below the waist. Here we discuss a mathematical model proposed to describe the growth of corneal ulcers caused by HSV Type I . This infection is common and frequently causes corneal opacification. Traditionally the morphology of the corneal ulcers has been described as either dendritic or amoeboid. The dendritic ulcers are by far the more frequent form and, though they are self-limited in general, ocasionally they can enlarge progressively changing to the amoeboid form. This is actually the natural course of the infection in the case of immunocompromised hosts or of inappropriate use of topic corticosteroids. In general, the amoeboid ulcers have a prolonged clinical course when compared to the dendritic ones. Regardless of their morphology, the ulcers are epithelial lesions that extend through the basement membrane whose swollen epithelial borders contain active viral particles.

In order to carry out a more quantitative study of the ulcer morphology, the fractal dimension of clinically diagnosed HSV ulcers (including both dendritic and amoeboid forms) have been estimated suggesting that their outlines are fractal objects . While the dendritic ulcers include branching and linear lesions, the geographic ulcers are no longer linear and as they increase in size, their perimeters become less and less irregular. In addition to its usefulness as a classification tool, the fractal properties of the ulcers may give information on the underlying mechanisms of viral spread within the epithelial tissue. For instance, a theory based only on the neurotropism of HSV and the dendritic-like distribution of nerve terminals can explain the branching pattern observed in dendritic ulcers , but it fails to explain the decrease of the fractal dimension (perimeter irregularity) with increasing ulcer sizes. An alternative explanation put forward by Landini et al , which will be the main focus of this paper, considers the ulcer shape as the natural outcome of the contiguous spread of viral particles modulated by variations in the cell susceptibilities to infection. To take into account the fact that viruses only infect cells that have appropriate receptor molecules on their surface, those authors proposed a cellular automaton model for the HSV I spread in which the corneal epithelial tissue is modeled by a two-dimensional lattice. In their model, each lattice site may be occupied either by a permissive cell (with probability $`q`$) or by a resistant cell (with probability $`1q`$). More pointedly, a permissive cell becomes infected whenever there are at least one infected cell in its neighborhood, while a resistant cell becomes infected if the number of infected and dead neighboring cells is larger than or equal to the integer parameter $`R>1`$ that measures the degree of resistance of the cell .

The simulated ulcers obtained with the cellular automaton have the same qualitative features of the clinical lesions and, in addition, for appropriate choices of the degree of resistance $`R`$ a dramatic change on the morphology of the ulcers is observed as the initial concentration of permissive cells $`q`$ increases beyond a certain value . This phenomenon was conjectured to be of a (qualitatively) similar nature as the ordinary percolation phase transition. The main contribution of this paper is to show, through the calculation of the dynamic and static critical exponents, that in the cases where a phase-transition does occur ($`R5`$), the transition belongs indeed to the universality class of the ordinary percolation . To carry out this analysis we use the so-called dynamic Monte Carlo method or spreading analysis  whose idea is to study the spreading of the infection starting from a configuration with a single infected cell on the center of the lattice. Clearly this technique is very well suited to our investigation since the characterization of the spreading behavior of the infection is exactly the issue we address in this paper.

The remainder of the paper is organized as follows. Following Landini et al , in Sec. 2 we give the set of rules that govern the evolution of the HSV I infection in a two-dimensional square lattice and present the evidences for the existence of a threshold phenomenon or phase transition for $`R5`$. In Sec. 3 we characterize this phase transition using the dynamic Monte Carlo method which allows the computation of the critical dynamic exponents that describe quantitatively the spreading of the infection from a single infected cell. Finally, some concluding remarks are presented in Sec. 4.

## 2 Model

The cellular automaton model is defined in a square lattice consisting of $`\left(L+1\right)\times \left(L+1\right)`$ sites, where each site is associated to a cell. Each cell is modeled by a four-state automaton corresponding to the different states of this cell: healthy permissive, healthy resistant, infected and dead. Except for the central cell, the initial state of any cell in the lattice is set either as permissive or resistant with probabilities $`q`$ and $`1q`$, respectively, so that there are no dead cells at the outset. The infection spreads from the single central infected cell and the ulcer (i.e., the cluster of dead cells) grows according to the following deterministic rules :

* An infected cell dies in the next time step.
* A healthy permissive cell becomes infected if at least one of its neighboring cells is infected.
* A healthy resistant cell becomes infected if at least $`R>1`$ of its neighboring cells are infected or dead.

The neighborhood of a given cell consists of its first and second nearest neighbors (Moore neighborhood). The infection and subsequent death of a resistant cell surrounded by $`R`$ or more dead cells is justified by the lack of tissue support. In addition, this is necessary to prevent the occurrence of large ulcers with small islands of resistant cells, which are not observed clinically . The four-state automaton considered allows transitions of one of the healthy states to the infected and dead states in a cyclic manner. At each time step we perform a parallel updating of all cell states.

For $`R8`$ we are dealing with a variant of the so-called diffusion percolation process where the geometry changes via a dynamic process and the nature of the growth depends on the local environment . For finite lattice sizes and open boundary conditions the above rules are repeated until either there are no more cells to infect or an infected cell reaches the lattice boundary. These different modes of termination generate dendritic (self-limited) and amoeboid (irrestricted) ulcers, respectively. It is interesting to note that the ordinary site percolation process is recovered for $`R>8`$, since in this case a resistant cell can never become infected and so the infection can propagate only through the permissive cells.

To illustrate the dependence of the different termination modes, and hence ulcer forms, on the control parameters $`R`$ and $`q`$ of the model we present in Fig. 1 the fraction of irrestricted ulcers generated in $`1000`$ runs. Each run corresponds to a different initial configuration of the lattice. Clearly, this fraction can be identified with the percolation probability $`\mathrm{\Pi }`$ . Rather interestingly, we have found that the results for $`R6`$ are indistinguishable within the numerical precision. Actually, this is expected since there is a preferred direction for the propagation of the infection, namely, from the center to the lattice boundaries, and so only the neighborhood facing the infection front matters to update the state of a healthy cell. Since the largest size of the error bars in this as well as in the next figure is twice the size of the symbols, they were ommited for the sake of clarity. More importantly, we have found that for $`R4`$, the results become independent of the lattice size already for $`L101`$.

However, for $`R5`$ the dependence on the lattice size, illustrated in Fig. 2 for $`R=5`$, indicates the occurrence of an atypical threshold phenomenon at a critical value $`q_c`$ in the limit $`L\mathrm{}`$. In fact, as $`q`$ increases from $`0`$ to $`1`$ the percolation probability $`\mathrm{\Pi }`$ vanishes for $`q<q_c`$, undergoes a discontinuous transition to some value $`\mathrm{\Pi }=\mathrm{\Pi }_c>0`$ at $`q=q_c`$ and then increases monotonically towards $`1`$. This transition is atypical in the sense that $`\mathrm{\Pi }_c`$ is not equal to $`1`$ above $`q_c`$, as in the case of the ordinary percolation transition , which means that in this regime there is a finite probability that the infection does not percolate, i.e., a dendritic ulcer is formed. The reason for that is due simply to the fact that the spreading process starts from a single central cell so that if the infection happens to percolate in a lattice of a given size then it is certain to percolate in a smaller lattice too, i.e., $`\mathrm{\Pi }\left(L_1\right)\mathrm{\Pi }\left(L_2\right)`$ for $`L_1<L_2`$. In particular, $`\mathrm{\Pi }\left(3\right)=1\left(1p\right)^8`$ yields an upper bound to $`\mathrm{\Pi }\left(\mathrm{}\right)`$. Of course, if the initial setting is such that there is an extensive number of infected cells, say $`\alpha L`$ with $`\alpha <1`$, randomly distributed over the bottom side of the lattice and periodic boundary conditions on the lateral sides, then the usual result $`\mathrm{\Pi }_c=1`$ is recovered . In fact, since the curves for different lattice sizes do not cross, the standard finite size scaling analysis aiming at determining both $`q_c`$ and the spatial correlation length exponent $`\nu _{}`$ for $`R5`$ (see, e.g., Ref. ) fails spectacularly and so we have to resort to other means to estimate those quantities.

## 3 Spreading analysis

We turn now to the analysis of the spreading behavior of the viral infection starting from a single infected cell located in the center of a lattice of infinite size. Finite size effects are absent because the lattice size is taken large enough so that during the time we follow the evolution the infection front can never reach the lattice boundaries. This of course sets an upper limit to the time we can follow the viral spread and so, for instance, for lattices of size $`L=4005`$ we let the infection evolve up to $`t=2000`$. As usual, we concentrate on the time dependence of the following key quantities : (i) the average number of dead and infected cells $`n\left(t\right)`$; (ii) the survival probability of the infection $`p\left(t\right)`$; and (iii) the average mean-square distance over which the ulcer has spread $`r^2\left(t\right)`$. For each time $`t`$ we carry out $`10^4`$ independent runs, hence $`p\left(t\right)`$ is simply the fraction of runs for which there is at least one infected cell in the lattice at time $`t`$. At the transition point $`q_c`$ we expect that the measured quantities obey the following scaling laws

$`p\left(t\right)`$ $``$ $`t^\delta `$ (1)
$`n\left(t\right)`$ $``$ $`t^\eta `$ (2)
$`r^2\left(t\right)`$ $``$ $`t^z`$ (3)
where $`\delta `$, $`\eta `$ and $`z`$ are dynamic exponents. Since the fractal dimension $`d_f`$ of the ulcer at a given time $`t`$ is defined as $`n\left(t\right)r^{d_f}`$ we have

$$d_f=2\frac{\eta }{z}$$
(4)
at the critical point. Note that this equation is different from the one used in the studies of directed percolation (see, e.g., Ref. ) because in the present case all runs generate an ulcer and so $`n\left(t\right)`$ as well as $`r^2\left(t\right)`$ are averages taken over all runs.

In Figs. 3, 4 and 5 we present log-log plots of $`p\left(t\right)`$, $`n\left(t\right)`$ and $`r^2\left(t\right)`$, respectively, as functions of $`t`$ in the vicinity of the critical point for $`R=5`$. The asymptotic straight lines observed in these figures are the signature of critical behavior while upward and downward deviations indicate supercritical ($`q>q_c`$) and subcritical ($`q<q_c`$) behaviors, respectively. We recall that in the subcritical regime only dendritic ulcers are formed, while in the supercritical regime the formation of amoeboid ulcers is much more frequent (see Fig. 2). The data shown in Fig. 3 yield $`q_c=0.3945\pm 0.0002`$ where the error is estimated by determining two values of $`q`$ as close as possible to the critical point for which upward and downward deviations can be observed. A precise estimate for the dynamic critical exponents is obtained by considering the local slopes of the curves shown in the previous figures. For instance, the local slope $`\delta \left(t\right)`$ is defined by

$$\delta \left(t\right)=\frac{\mathrm{ln}\left[p\left(t\right)/p\left(t/8\right)\right]}{\mathrm{ln}8},$$
(5)
which for large $`t`$ behaves as

$$\delta \left(t\right)\delta +\frac{a}{t}$$
(6)
where $`a`$ is a constant. Analogous expressions hold for $`\eta \left(t\right)`$ and $`z\left(t\right)`$. Hence plots of the local slopes as functions of $`1/t`$ allow the calculation of the critical exponents. Applying this procedure for the critical curves we find the exponents $`\delta =0.0870\pm 0.0001`$, $`\eta =1.5866\pm 0.0007`$ and $`z=1.6843\pm 0.0003`$. The errors in the critical exponents are, as usual, the statistical errors obtained by fitting the local slopes by straight lines in the large $`t`$ regime. We expect, of course, that the (uncontrolled) systematic errors are much larger than those. Using Eq. (4) we obtain $`d_f=1.8840\pm 0.0005`$ which is in very good agreement with the analytical prediction for the ordinary percolation $`d_f=91/481.896`$ .

As the dynamic exponents $`\delta `$, $`\eta `$ and $`z`$ for the ordinary percolation problem are not very well known, to show unambiguously that this ulcer formation model belongs to the universality class of the ordinary percolation we ought to estimate the static exponents $`\beta `$ and $`\nu _{}`$. We recall that the exponent $`\beta `$ gives a measure of how the fraction of lattice cells belonging to an infinite cluster vanishes as the percolation threshold $`q_c`$ is approached in the supercritical regime while $`\nu _{}`$ is the correlation-length exponent in the space direction. To do so we calculate first the exponent $`\nu _{}`$ which governs the decay of the concentration of infected cells $`i\left(t\right)`$ in the subcritical regime. In fact, since in this regime the correlations are short-ranged one expects $`i\left(t\right)`$ to decay exponentially

$$i\left(t\right)A\left(q\right)\mathrm{exp}\left[\left(q_cq\right)^\nu _{}t\right]t\mathrm{}$$
(7)
where $`A\left(q\right)`$ is some time independent function. Fig. 6 not only illustrates the adequacy of this assumption but permits also the evaluation of the decay constant

$$\lambda =\left(q_cq\right)^\nu _{}$$
(8)
from the asymptotic slopes of the curves $`\mathrm{ln}i`$ vs. $`t`$. The results presented in Fig. 7, showing the dependence of $`\lambda `$ on the distance $`q_cq`$ from the critical point, allows the calculation of the exponent $`\nu _{}`$ as the slope of the straight line, yielding $`\nu _{}=1.54\pm 0.03`$. Once this exponent is known we can use the scaling relations $`\beta =\nu _{}\delta `$ and $`\nu _{}=z\nu _{}/2`$  to estimate the static exponents. We find $`\beta =0.134\pm 0.003`$ and $`\nu _{}=1.30\pm 0.03`$ which, within error bars, are in agreement with the exact values of the corresponding exponents of the ordinary percolation, namely, $`\beta =5/360.139`$ and $`\nu _{}=4/31.333`$ .

We have carried out a similar spreading analysis for $`R6`$ and, as hinted in Fig. 1, we have found a slightly larger percolation threshold, namely, $`q_c=0.4075\pm 0.0002`$ which, within error bars, is shown to be independent of the value of $`R6`$. Furthermore, since the larger the resistance parameter $`R`$, the more similar the ulcer formation problem is to the ordinary site percolation, we have found the same dynamic and static critical exponents as for the case $`R=5`$, as expected.

## 4 Conclusion

Using the dynamic Monte Carlo method we have shown unambiguously that the phase transition observed in the model for formation of herpes simplex ulcers proposed by Landini et al  belongs to the universality class of the ordinary percolation. The value of this finding should not be underrated since the infection process actually resembles a diffusion percolation process where the growth depends on the local environment, in the sense that the decision on whether or not a resistant cell will become infected depends on the time-dependent states of several of its neighbors. Furthermore, since the ulcer formation model described here may be thought of as a damage spreading process, one could expect that the transition were in the universality class of the $`\left(2+1\right)`$ directed percolation instead. However, as pointed out by Grassberger , this is not so because in the ulcer formation model the damage never heals (even if it does not spread), i.e., the probability that an infected or dead cell becomes healthy is zero.

The finding that for $`R<5`$ the model does not present a phase transition reflects the nontrivial role played by the resistance parameter $`R`$ in this percolation process. In these noncritical cases the probability that an infinite or irrestricted ulcer is generated is given by the smooth size-independent curves shown in Fig. 1. Similarly to a noncritical forest fire model  the growing of this type of ulcer may be characterized by infection fronts with fractal dimension $`D`$ whose value probably depends on the resistance parameter $`R`$ (of course, $`D=1`$ for $`R=1`$). An additional feature that makes the quantitative study of this viral spreading model rather challenging is the result that the percolation probability curves for different lattice sizes do not cross (see Fig. 2), which complicates enormously the estimate of the percolation threshold and critical exponents through the standard finite size scaling method.

Some remarks on the biological interpretation of our results are in order. According to the specialized literature , amoeboid ulcers are, in general, observed in immunocompromised patients or in patients that made inappropriated use of corticosteroids. In the present model, these conditions would correspond to a decrease of the degree of resistance $`R`$ of the resistant cells or to an increase of the initial concentration $`q`$ of permisssive cells. Although this model does not take into account the recurrent characteristic of this kind of infection, in which case the variability of $`q`$ would probably play an important role, nor the possibility of variation of $`R`$ during the course of the infection, its predictions are in qualitative agreement with the clinical observations. In fact, Fig. 1 points out the prevalence of amoeboid ulcers when $`R`$ decreases or $`q`$ increases. This agreement lends support to the hypothesis that the morphology of the ulcers is determined by the viral spreading through cells with different susceptibilities to infection.

To conclude we should mention that an extension of the original model proposed by Landini et al in which both the regeneration of dead cells as well as the spontaneous outbreak of infection anywhere in the lattice are taken into account has already been considered in the literature . Interestingly, in this case the viral spreading model becomes very similar to the critical forest fire model with immune trees . In particular, the resistance parameter $`R`$ of the ulcer formation model is akin to the immunity probability, i.e., the probability that a tree is not ignited though one of its neighbors is burning. According to a conjecture put forward by Grassberger , the extended ulcer formation model should be in the universality class of directed percolation, since it allows for the regeneration of dead cells. This suggestion is strengthened by the finding that the forest fire model with immune trees is in that class of universality .

The work of J.F.F. is supported in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Proj. No. 99/09644-9. The work of R.M.Z.S. is partially supported by CNPq. C.P.F. is supported by FAPESP. We thank FAPESP for supporting R.M.Z.S’s visit to São Carlos where part of her work was done.


text:
# Second-parameter Globular Clusters in the Milky Way and in M33 as Tracers of Mass Loss from M31 in the Early Epoch?

## 1. Problem and Main Questions

There are contradictory conclusions about the dwarf spheroidal (dSph) galaxies’ tidal dissolution and its rate, as well as a number of observational constraints on this process in the Galaxy (e.g., van den Bergh 1996; Gilmore et al. 1999; Shetrone et al. 2000). Moreover, preliminary results on the CMD study of GCs in massive spirals of LG suggest a fraction of SPGCs to be at least low among GCs in M31 (Rich et al. 2001) and to be high as for M33 GC population (Sarajedini et al. 1998). Hence there appear doubts about tidal disruption/stripping of dSphs as universal mechanism of and main contributor to a formation of populations of the SPGCs in the above three galaxies. Note also that Chandar et al. (2001) have estimated the total number of GCs in M 33 to be $`75\pm 14`$ which gives a specific frequency significantly higher than in other late-type spirals.

By now, it is established that the Galactic population of SPGCs shows a number of basic distinctive characteristics as compared with those of other GCs in the halo of the Milky Way (e.g., Da Costa 1994; van den Bergh 1996). Some of these characteristics imply that the Galactic SPGCs are very likely to have an accretive nature. However, if a tidal dissolution of dSphs is not main and the only contributor to the formation of population(s) of the SPGCs in our Galaxy (and in M33) then the following vexed questions arise.

What is a material accreted and where did it originate and come from? Why has the formation of SPGCs lasted after formation of the bulk of halo objects in rapid collapse?

## 2. Putative Role of the Mass Outflows

The high-redshift Layman break galaxies (LBGs) are suggested to be progenitors of the present-day massive spheroids being in particular components of luminous early-type spirals (Friaca & Terlevich 1999). As summarized by Heckman (2000) starburst-driven galactic superwinds observed in LBGs play important role in the mass outflow from them. During main, most powerful episode ($`\mathrm{𝟐}`$ Gyr) of formation of the Pop. II stars in a star forming galaxy like M31, its outflows may carry off the galaxy the total mass of gas as large as mass of the stellar bulge and spheroid, with typical speeds ranging from a few $`\mathrm{𝟏𝟎}^\mathrm{𝟐}`$ to $`\mathrm{𝟏𝟎}^\mathrm{𝟑}`$ km/s (Pettini et al. 2001). Some a few tens percent of this mass ($`𝚫𝐌\mathrm{𝟏𝟎}^{\mathrm{𝟏𝟏}}𝐌_{}`$ or so) with velocities exceeding the escape velocity of M31 could have quited the galaxy, portion of which might be accreted on our Galaxy (and M33). If we accept present distance between the Galaxy and M31 ($`\mathbf{0.7}`$ Mpc), ratio of their mass ($`𝐌_𝐆/𝐌_𝐀\mathrm{𝟏}/\mathbf{1.5}`$), and also a mean velocity of expanding gas to be of order $`\mathrm{𝟒𝟎𝟎}`$ km/s, then we easily find that a time spent by gas to reach a region of the Lagrangian point between the two galaxies, may be as long as $`\mathrm{𝟏}`$ Gyr. An upper limit of mass $`\delta 𝐌`$ of gas accreted on the Galaxy is estimated from a formula obtained (in the approach of gas accretion in gravity potential $`\phi R^1`$ of a field of less massive companion in a binary system) and kindly offered by Postnov (2000, private communication):

$$\delta M\frac{1}{4}\left(\frac{M_G}{M_A}\right)^2\mathrm{\Delta }M\frac{1}{10}\mathrm{\Delta }M10^{10}M_{}$$

Even if a mass of really accreted gas was an order of magnitude lower than the upper limit estimated and if in turn approximately ten percent of this mass was converted into objects, i.e., $`\mathrm{𝟏𝟎}^\mathrm{𝟖}𝐌_{}`$, that would be equal to a total mass of the bulk of SPGCs, of a few dSphs, and of some portion of field stars which may have formed due to such a process in the outer halo of our Galaxy.

Similarly, it is possible to conclude that objects assumed to form from a gas accreted onto M33 might comprise a total mass of order a few $`\mathrm{𝟏𝟎}^\mathrm{𝟔}𝐌_{}`$. Some 10 to 20 GCs (with fraction of field stars formed) with cluster’s typical mass of order $`10^5M_{}`$ may amount just to this mass.

### Acknowledgments.

I am deeply thankful for the IAU travel grant. Also, I am very grateful to Dr. A.A. Suchkov for his encouraging and valuable comments on the scenario presented, and to Drs. K. Postnov, Yu.A. Shchekinov, Yu.N. Efremov, N.I. Shakura, and M.V. Sazhin for helpful discussions.

## References

Chandar, R., Bianchi, L., & Ford, H. C. 2001, A&A, 366, 498

Da Costa, G. S. 1994, In ESO Conference and Workshop Proc., 51, The Local Group: Comparative and Global Properties, ed. A. Layden, R.C. Smith & J. Strom, 101

Friaca, A. C. S. & Terlevich, R. J. 1999, MNRAS, 305, 90

Gilmore, G., Hernadez, X., Valls-Gabaud, D. 1999, astro-ph/9910409

Heckman, T. M. 2000, astro-ph/0009075

Pettini, M., Steidel, C. C., Shapley, A. E., et al. 2001, astro-ph/0101254

Rich, R. M., Corsi, C. E., Cacciari, C. et al. 2001, this volume

Sarajedini, A., Geisler, D., Harding, P., & Schommer, R. 1998, ApJ, 508, L37

Shetrone, M., Cote, P., Sargent, W. L. W. 2000, astro-ph/0009505

van den Bergh, S. 1996, PASP, 108, 986