format: text

text:
# GZ Cancri: a cataclysmic variable at the lower edge of the period gap Based on observations made at the European Southern Observatory

## 1 Introduction

Cataclysmic variables (CVs) are close interacting binaries with a white dwarf as primary component and usually a late-type main-sequence star as secondary. In the course of their evolution from wide detached binaries, angular momentum loss thought to be driven by magnetic braking and gravitational radiation has shrunk the Roche lobe of the secondary into contact with the stellar surface. This causes mass transfer via the Lagrange point $`L_1`$ into the gravitational well of the primary. In the absence of strong magnetic fields, an accretion disc around the white dwarf is formed, which is also the major source of the typical emission line spectrum. For a comprehensive overview on CVs see Warner (warn95 (1995)).

GZ Cnc (= RX J0915.8+0900 = Tmz V034) was discovered independently as a variable star and a ROSAT source (Takamizawa taka98 (1998); Bade et al. bade+98 (1998)). A spectrum was presented by Jiang et al. (jian+00 (2000)), showing a blue continuum and moderately strong emission lines. Although the spectral coverage went up to 8200 Å, no spectroscopic signatures of the secondary star were detected. Kato et al. (kato+01 (2001)) obtained lightcurves of the star during decline in outburst. Their longest time series spans $``$3.1 h and does not show any apparent periodic variations. Both the spectroscopic and the photometric properties, together with the outburst behaviour (slow rise, small amplitude; Kato et al. kato+01 (2001)) therefore point to an SS Cyg-type dwarf nova with an orbital period $`>`$ 3 h. In a recent paper Kato et al. (kato+02 (2002)) investigate the long-term lightcurve in more detail, and emphasise similarities to those of intermediate polars. Again, this would favour a long orbital period, as the vast majority of those systems is found above the period gap (Ritter & Kolb rittkolb98 (1998)).

We here present time-resolved spectroscopy in order to derive the orbital period, which, in fact, turns out to be less than 3 h. In Sect. 2 we describe the observations and the reduction process. The results are presented in Sect. 3 and subsequently discussed in Sect. 4. We complete our paper with a short conclusion in Sect. 5.

## 2 Observations and reductions

The data were taken on January 16–18, 2001, and on March 26–28, 2002, with the DFOSC system at the 1.54 m Danish telescope at ESO, La Silla, Chile. The spectroscopy was performed using grism #7, yielding a wavelength range of 3800–7000 Å. A slit width of 1$`\stackrel{}{.}`$5 gave a spectral resolution of 4.9 Å. The measurements were accompanied by He-Ne arc spectra for wavelength calibration and by spectrophotometric standards for flux calibration. Between each row of spectroscopic measurements, photometric exposures were taken in the $`V`$ passband. Additionally, calibrated $`BVR`$ photometry was obtained on 2002-03-26. Table 1 presents the observation log.

The reduction process included bias and flatfield correction, as well as wavelength and flux calibration, which were performed with IRAF routines (Tody tody93 (1993)). Further analysis of the spectroscopic data was done with ESO-MIDAS (Warmels warm92 (1992)), while the photometry made use of the IRAF daophot package and the standalone daomatch and daomaster routines (Stetson stet92 (1992)). The photometric calibration was obtained by comparison with Landolt (land92 (1992)).

## 3 Results

### 3.1 Photometry

Only on 2002-03-26 direct calibrated photometric data were obtained in $`B`$, $`V`$, and $`R`$ passbands. The resulting colours are

$`BV=0.097(99),`$ $`VR=0.394(97),`$
$`\mathrm{at}V=15.471(85),`$
which lie in the usual range occupied by dwarf novae (Echevarría eche84 (1984)), and put it clearly away from the main sequence (Fig. 1).

For the other nights, calibrated $`V`$ magnitudes for GZ Cnc were computed with respect to an average lightcurve composed of 5 stars which could be measured in all frames and proved to be constant within the errors. The resulting magnitudes are given in Table 2. During the 2002 observing run the system appears to have been in a stable quiescent state, with variations $`\mathrm{\Delta }m_V0.3`$ mag probably due to flickering or an additional continuum source (e.g., the bright spot). A comparison with Kato et al. (kato+02 (2002)) shows that our observations lie close to two recorded outbursts on 2002-03-11 and 2002-04-02. The 2001 data, however, clearly show a brightness increase, which could be interpreted as the onset of an outburst (see Fig. 2), although to our knowledge none has been registered in the subsequent days. However, the coverage of GZ Cnc contains a lot of gaps. For the year 2001, only 11 observations are recorded in the VSNET database111The archive of the VSNET mailing list is located at http://www.kusastro.kyoto-u.ac.jp/vsnet/Mail, among them one single positive detection. A weak outburst could thus easily have been missed.

The increase rate of roughly 0.2 mag d$`{}_{}{}^{}1`$ in our data fits well with the outburst recorded by Kato et al. (kato+01 (2001)), that showed a rate of $``$0.25 mag d$`{}_{}{}^{}1`$ over 5 nights. This lends further support to the idea of an overlooked outburst in 2001.

### 3.2 Properties of the average spectra

The spectrum of GZ Cnc (Fig. 3) resembles that of a dwarf nova, with dominating emission lines from the Balmer and He I series. Further lines include He II $`\lambda `$4686, as well as Fe II $`\lambda `$5172. Table 3 summarises the identified lines and their equivalent widths. The emission lines are of rather moderate strength (e.g., Williams will83 (1983)), which indicates a moderately high mass-transfer rate (Fig. 6 in Patterson patt84 (1984)). As noted by Kato et al. (kato+02 (2002)), He II appears to be somewhat stronger than in most dwarf novae, and especially those with a short orbital period.

The slow rise in the photometry during the 2001 run is also roughly reflected by the equivalent widths of the emission line, which first show a slight rise from the first to the second night, followed by a steep decline on the third night (Table 3, Fig. 3). It appears thus that this variation is indeed due to a brightening of the disc continuum.

### 3.3 Period analysis

The Doppler variation of the emission lines caused by the orbital motion of the system has been measured with a Gaussian fit to the whole H$`\alpha `$ line profile. Radial velocities were computed with respect to the O I night sky emission line at $`\lambda `$5577 Å. As the mean velocity showed strong nightly variations, and the respective time ranges proved to cover more than one orbit, we normalised the velocity data individually by subtracting the respective mean velocities.

The search for periodicities was conducted using the Scargle (scar82 (1982)) and AOV (Schwarzenberg-Czerny schw89 (1989)) algorithms as implemented in MIDAS’ tsa context. Unfortunately, the time span between both observing runs proved to be too long for an analysis of the combined data, so that the runs had to be investigated individually. The Scargle algorithm produced five significant peaks in the range $`f=9.3513.35\mathrm{cyc}\mathrm{d}^1`$. The frequencies corresponding to the highest peaks were 10.34 for the 2001 data and 11.34 for 2002. The AOV method in both cases picked out $`f11.35`$. For the 2002 data, this peak is well separated from the next highest one, for 2001 the two peaks of 10.34 and 11.34 are very close (Fig. 4, lower and upper plot).

In most CVs, additional emission components distort the symmetric, disc-borne, line profile (Tappert & Hanuschik tapphanu01 (2001)), with the type and strength of this distortion being likely to depend on the luminosity state of the system (e.g., Tappert tapp99 (1999)). Since on 2001-01-18 both the photometric magnitudes (Table 2) and the emission line equivalent widths (Table 3) differ strongly from those of the other nights in this run, we also conducted a period analysis for the data of the first two nights in 2001 only. Indeed the peak at $`f=11.3`$ in the AOV periodogram is now much more pronounced with respect to the other frequencies (Fig. 4, middle plot).

In order to test this period against the two most prominent aliases at 10.3 and 12.3 cyc d$`{}_{}{}^{}1`$ (referring to the 2002-03 data) we have run a Monte Carlo simulation (Mennickent & Tappert menntapp01 (2001)). This routine computes a thousand data sets by adding a random value to each data point, based on the variance of the sine fit to the real data. A periodogram for each of these sets is established and the frequencies and their power are recorded. Equivalently to the method outlined by Thorstensen & Freed (thorfree85 (1985); hereafter TF85) one can now compute a discriminatory power $`p_d`$ by dividing the number of times where a specific frequency has the largest power in comparison to the other aliases by the number of total tries. The second statistical measure used by TF85, the correctness likelihood $`p_c`$, can be obtained by counting the number of times when the ratio of the power of the highest peak and that of the second highest one for the artificial data corresponds to that of the real data within a certain range of ratios (TF85 propose a typical value of 0.1).

Table 4 summarises the results of this analysis for the three most prominent frequencies. The highest statistical significance for the principal frequency $`f_211.3`$ cyc d-1 is obtained to $`p_c`$ = 0.65 in the 2002-03 data. This favours the identification of $`f_2`$ with the orbital frequency, but is certainly too low to provide an unambiguous decision, with $`f_110.3`$ cyc d-1 representing an important alias.

After the initial submission of this paper we have been informed of the research by Fenton & Thorstensen (private communication) who obtained $`P=0.0880(2)`$ d from time-resolved spectroscopy of GZ Cnc. They report that their run spans over 6 h of hour angle and does not present any daily alias problems.

This period fits very well with our best frequency $`f_2`$, which, as an average of the ‘quiescence’ data sets, corresponds to a period

$$P=0.08825(28)\mathrm{d}=2.118(07)\mathrm{h}.$$

We therefore identify this period with the orbital variation.

Fig. 5 presents the radial velocities folded on our period.

### 3.4 Line-profile analysis

Even for medium spectral resolution an analysis of the emission line profile is always worth undertaking, as at least it might provide a motivation for further studies. In the present case, the fact that at least one data set caught the system in a different luminosity state, and that its exclusion led to a more unambiguous determination of the orbital period, already calls for a more thorough investigation. The only really feasible line in our data for such a study is H$`\alpha `$, as other lines are either too noisy (He I, He II, etc.) or are distorted by blends (basically all Balmer lines from H$`\beta `$ on bluewards).

Due to the modest spectral and phase resolution the more advanced technique of Doppler tomography (Marsh & Horne marshorn88 (1988)), which in principle would yield the clearest picture on the emission distribution in the system, does not provide any significant results. However, the application of the diagnostic diagram (Shafter shaf83 (1983)) and the $`V/R`$ plot (Tappert tapp99 (1999)) will still give at least a rough idea of a possible additional component and of the differences between the individual data sets.

In the former method, two identical Gaussians (in the present case with FWHM = 4 Å), horizontally separated by a quantity $`d`$, are used to measure radial velocities of different parts of the emission line by variation of $`d`$ (Schneider & Young schnyoun80 (1980)). The resulting velocities $`v(d,\phi )`$ are fit with a function

$$v(d,\phi )=\gamma (d)K_1(d)\mathrm{sin}(\phi ).$$

The fit parameters are plotted as a function of $`d`$ to yield the diagnostic diagram. In order to identify the part of the line which is undisturbed from possible additional components, one usually also plots the percentage error of the semi-amplitude $`K_1`$, $`\sigma (K_1)/K_1`$. A sharp increase in the latter at higher separations can be interpreted as the continuum noise dominating over the signal from the emission line. The last separation before this happens, $`d_{\mathrm{max}}`$ is usually taken as that part of the line that reflects best the motion of the white dwarf, especially if the corresponding parameters tend to approach constant values. Note, that this assumes that the wings of the line profile stem from material which is symmetrically distributed within the accretion disc.

With the comparatively low resolution of our data we do not expect to be able to extract significant parameters from the diagnostic diagram, but are mainly interested in comparing the behaviour of the individual data sets. We have first computed diagrams for each night separately. In the case of the sets from 2001-01-16 and 2001-01-17, and in the case of 2002-03-27 and 2002-03-28, the results were identical within the errors, and we therefore used the combined data sets. The resulting diagrams are shown in Fig. 6. The $`\sigma (K_1)/K_1`$ curves do not really give a good indication for $`d_{\mathrm{max}}`$, and thus our choice to correct the zero point of the phase for the value corresponding to $`d`$ = 60 Å is somewhat arbitrary. We were originally motivated by the behaviour of the 2002-03 $`K_1`$ and $`\sigma (K_1)/K_1`$ curves, and for the 2001-01-16+17 data by the fact that at this point the $`K_1`$ value corresponds well to the 2002-03 $`K_1`$. However, we are not particularly interested in a specific value for $`\gamma `$ or $`K_1`$, and although $`\phi _0`$ is important for the interpretation of the subsequent $`V/R`$ plot, it stays almost constant over the whole range of separations, and its choice therefore does not appear to be critical.

From the appearance of the curves in the diagnostic diagram we can derive two pieces of information. First, that the variation of $`K_1`$ provides evidence for the presence of an additional emission component, and second, that the similarity of the variation in all parameters of the 2001-01-16+17 and 2002-03 data sets indicates a similar type of additional emission, while there are strong differences with respect to the 2001-01-18 data.

For the $`V/R`$ plot the fluxes of the blue (violet) half, $`F(V)`$, and of the red half, $`F(R)`$, of the emission line profile are computed to obtain

$$V/R=\mathrm{log}\frac{F(V)}{F(R)}.$$

The outer borders of the line profile are defined by a certain intensity level $`I`$, the central wavelength which separates both halves is taken as the mean of the two outer wavelengths $`\lambda _V(I)`$ and $`\lambda _R(I)`$. For the present case we chose $`I=0.15I_{\mathrm{max}}`$, with $`I_{\mathrm{max}}`$ being the peak value of the emission line.

The resulting plots, $`V/R`$ vs. phase, are given in Fig. 7. While the variation is very noisy and badly defined in the 2001-01-16 and 17, and the 2002-03 data, they do show a similar behaviour. In both cases a maximum at similar phases (0.9 in the upper plot, 1.0 in the lower one) can be identified. As the respective zero points of the phases were determined individually, we do not consider the observed phase difference as evidence for a fundamentally different line profile variation. The plot for the 2001-01-18 data (middle plot), however, shows a much stronger, better defined, variation, which additionally is significantly shifted in phase.

Maximum and minimum in the $`V/R`$ plot can usually be identified as the maximum blue- and redshift of an additional emission component. This assumes that this component is dominating the line profile and mostly confined to the centre of the profile. However, this does not seem to work for GZ Cnc. Fig. 8 shows that the minimum in the 2001-01-18 data around phase 0.65 is due to a small bump in the blue wing which then moves redwards with phase. Note that the stronger part of the line profile shows very similar velocities in both the 2001-01-18 and the 2001-01-17 data.

While this additional component in H$`\alpha `$ at first glance appears to be an emission feature, a closer investigation of all lines (Fig. 9) and a comparison with respect to unaffected phases (Fig. 10) reveals that it is actually a high velocity absorption component. It is visible in all emission lines from phase 0.54 to 0.93, with the probable exception of He II.

As the component is most clearly seen in the He I $`\lambda `$5876 line, we chose the latter to measure its velocities. Fig. 11 shows that it appears at phase 0.54 at a velocity $`650`$ km s-1, moving redwards and apparently reaching a plateau around phase 0.9. The latter would suggest a confinement to negative velocities, but this impression is mainly caused by the first and the last data point, which certainly are the most uncertain ones (see Fig. 9). On the other hand, we find no evidence that the component ever appears at positive velocities (Fig. 12). Unfortunately, our data does not cover a second orbit, so that it remains unclear if the absorption component represents a transient or a periodic event.

## 4 Discussion

Our analysis revealed a period of 2.118(7) h, placing GZ Cnc just at the lower edge of the period gap of CVs (Fig. 13). This is a region which is sparsely populated, especially by dwarf novae, and it is therefore of particular interest to determine the CV subtype of the system.

Kato et al. (kato+02 (2002); hereafter K2002) studied the long-term behaviour, and found a similarity to those of intermediate polars (IPs) in that the outburst distribution appeared to be ‘clustered’. Specifically the authors mention the detection of four outbursts within two months in 2002, while the only known outburst before that date has been recorded at the end of 2000. Still, our data proves that at least one brightening has occurred during 2001 that remained undetected. It is therefore possible that the observed clustering is not due to a variation of the outburst frequency, but caused by a variable outburst amplitude, especially as the upper detection limit in the observations used by K2002 varies between $`m_{\mathrm{vis}}`$14.8 and 13.2, the latter being already very close to the object’s maximum intensity (Fig. 1 in K2002). Note, however, that frequent outbursts of variable amplitude by no means contradict the suggested similarity to V426 Oph: it certainly fits well with the latter’s long-term lightcurve (Fig. 2 in K2002).

Another point in favour of the IP scenario might be provided by the appearance of the absorption component during rise to outburst. Such features during certain phases have been observed e.g. in early outburst of the short period ($`P_{\mathrm{orb}}`$ = 1.63 h) IP EX Hya (Hellier et al. hell+89 (1989)) and in the long period ($`P_{\mathrm{orb}}`$ = 4.85 h) IP FO Aqr (Hellier et al. hell+90 (1990)). They are furthermore a common property of the SW Sex stars (Hellier hell00 (2000), and references therein), which have been recently proposed to represent high mass-transfer IPs by Rodríguez-Gil et al. (rodr+01 (2001)).

These features have been explained by a scenario where part of the gas stream from the secondary overflows the disc, and impacts on the inner part of the disc or connects with the magnetosphere of the white dwarf (Hellier hell00 (2000); Hellier et al. hell+89 (1989)). However, an important difference of the GZ Cnc absorption with respect to the above systems is that their absorption features consist of low-velocity components. A stellar wind as suggested for a similar phenomenon in the nova-like BZ Cam (Patterson et al. patt+96 (1996)) might represent a more probable alternative.

Another similarity with EX Hya, but also a common feature of dwarf novae, lies in the presence of an isolated emission component in the quiescence data, as suggested by the $`V/R`$ plot (Fig. 7). This component is confined to the line centre and could thus stem from the impact region of the gas stream with the outer accretion disc.

## 5 Conclusions

There are three pieces of evidence which might favour GZ Cnc as an intermediate polar (IP): the long-term behaviour (Kato et al. kato+02 (2002)), the unusual, although not exorbitant, strength of the He II line, and the appearance of the absorption component during rise to outburst. None of these is really conclusive, basically because they are based on insufficient data. However, in their sum they represent at least a strong indication.

We conclude that further observations are necessary in order to decide pro or contra an IP nature for GZ Cnc. Especially helpful would be a better sampled long-term lightcurve and high-speed photometry. If an IP nature were to be confirmed, it would make GZ Cnc a rare system indeed: the TPP catalogue (Kube et al. kube+02 (2002)) records only one out of twenty confirmed IPs below the period gap.

Finally, the non-uniform emission distribution in GZ Cnc marks it as an interesting system for further and more detailed studies using higher resolved spectroscopic data.

###### Acknowledgements.

We thank Boris Gänsicke and Ronald Mennickent for advice and discussions. John Thorstensen kindly informed us on the result of Bill Fenton’s and his research, which is gratefully acknowledged. We thank the referee, Pablo Rodríguez-Gil, for helpful comments and for pointing out the behaviour of high-velocity components in intermediate polars and SW Sex stars. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.


text:
# Collapses and explosions in self-gravitating systems

## I Introduction

Systems of particles interacting via a potential with attractive nonintegrable large $`r`$ asymptotics, $`U(r)r^\alpha `$, $`0<\alpha <3`$, and sufficiently short-range small $`r`$ regularization exhibit gravitational phase transition between a relatively uniform high energy state and a low-energy state with a core-halo structure pr ; ki2 ; ch1 ; usg ; dv ; ch ; chs ; chi . Extensive mean-field (MF) studies of the equilibrium properties of such systems pr ; ki2 ; ch1 ; usg ; dv ; ch ; chs ; chi  revealed that in a microcanonical ensemble during such a transition entropy has to undergo a discontinuous jump from a state that just ceases to be a local entropy maximum to a state with the same energy but different temperature, which is the global entropy maximum. Due to the long-range nature of gravitational interaction, the MF studies are believed to provide asymptotically (in the infinite system limit) exact information about the density and the velocity distributions and other thermodynamical parameters of the uniform state. The applicability of the MF theory to the description of the core-halo state is less obvious as the properties of a core are controlled by the short-range asymptotics of the potential.

Relatively little is known about how such a transition actually occurs, however. Youngkins and Miller bm2  performed a Molecular Dynamics (MD) study of a one-dimensional system consisting of concentric spherical shells. Their main emphasis was to check the MF description of the stable and metastable states rather than to study the dynamics of the phase transition itself. Cerruti-Sola, Cipriani, and Pettini pet  studied the phase diagram of a more realistic 3-dimensional particle system by using Monte Carlo and MD methods. Their studies again focused on the equilibrium properties of the system rather than on the dynamics of the transitions. In addition, their general conclusions about the second order of the gravitational phase transition apparently contradict the MF results ki2 ; usg ; dv ; chi .

Here, we attempt to resolve this contradiction. A MF description of the dynamics of collapse in ensembles of self-gravitating Brownian particles with a bare $`1/r`$ interaction based on a Smoluchowski equation was developed by Chavanis et al chs . It predicts a self-similar evolution of the central part of density distribution to a finite-time singularity. However, the precise nature of the random force and friction terms in the corresponding Fokker-Plank equation as well as the applicability of the overdamped limit used to reduce the Fokker-Plank equation to the Smoluchowski equation are not entirely understood. A more rigorous approach based on the Fokker-Plank equation with Landau collision integral was used by Lancellotti and Kiessling ki3  to prove a scaling property of the central part of the density profile. The model considered there allows the particles to escape to infinity and therefore does not have an equilibrium or even a metastable state.

There exists a vast amount of literature on cosmologically- and astrophysically- motivated studies of the temporal evolution of naturally occurring self-gravitating systems (see, e.g., Ref. cos  and references therein). The selection of systems and their initial and final conditions made in such studies are typically astrophysically-motivated; the considered systems are often too complex to make a general conclusions about the phase diagrams and phase transitions in such systems.

In this paper we present MD studies of gravitational collapse, and reverse to collapse, i.e., explosion, transitions in a microcanonical ensemble of self-attracting particles. Besides their pure statistical mechanical implications, these studies represent our attempt to bridge a gap between the usually complicated MD and hydrodynamic simulations of the realistic astrophysical systems and the MF analysis of the phase diagram of simple self-gravitating models.

A system with soft Coulomb potential $`(r^2+r_0^2)^{1/2}`$ is considered. Such systems are have been studied using both MF theory (see, e.g., Refs. usg ; chi ) and simulations pet . We chose the microcanonical ensemble as the most fundamental one for the long-range interacting systems. It has to be noted that the considered system is strongly ensemble-dependent: While the nature of the uniform state is the same in both microcanonical and canonical ensembles (apart form the difference in their stability range), the core-halo states and the collapse itself in these ensembles have very little in common with each other dv ; chs .

A MF phase diagram of the considered self-attracting microcanonical system is presented in Fig. (1) usg ; chi .

High- and low-energy branches terminating at the energies $`ϵ_{coll}`$ and $`ϵ_{expl}`$ correspond to the uniform and core-halo states. The energy $`ϵ^{}`$ where the entropies of the core-halo and uniform states are equal is the energy of the true phase transition; the uniform and the core-halo states are metastable in the energy intervals $`(ϵ_{coll},ϵ^{})`$ and $`(ϵ^{},ϵ_{expl})`$, respectively. However, for the phase transition to occur at or near $`ϵ^{}`$, a macroscopic-scale fluctuation with prohibitively low entropy is required. Consequently, the metastable branches are stable everywhere except at the vicinity of $`\mathrm{\Delta }ϵN^{2/3}`$ of their end-points $`ϵ_{coll}`$ and $`ϵ_{expl}`$ ka2 ; chi , where $`N`$ is the number of particles in the system. Hence it is natural to assume that once the energy of the system in the uniform state is set sufficiently near above or below $`ϵ_{coll}`$, the system will undergo a collapse to a core-halo state with the same energy and higher entropy. Similarly, if the energy of the core-halo system is set sufficiently near below or above $`ϵ_{expl}`$, the system will undergo an explosion bringing it to a uniform state with the same energy and higher entropy. Our goal here is to study if and how such collapse and explosion proceed in a realistic three-dimensional N-particle dynamical system.

The paper is organized as follows. In the next section we formally introduce the system, outline the MF analysis and describe the MD setup. Then, we present the simulation results for the equilibrium uniform and core-halo states and compare them to the MF predictions. After that we describe and interpret the observed dynamics of the collapse and the explosion transitions. A discussion of the obtained results concludes the paper.

## II simulation

We consider a system consisting of $`N`$ identical particles of unit mass confined to a spherical container of radius $`R`$ with reflecting walls. The particles interact via the attractive soft Coulomb pair potential $`(r^2+r_0^2)^{1/2}`$. Using a traditional convention for self-gravitating systems in which the equilibrium properties of such systems become universal, we define rescaled energy $`ϵ`$, inverse temperature $`\beta `$, distance $`x`$, velocity $`u`$, and time $`\tau `$ as

$`ϵE{\displaystyle \frac{R}{N^2}}`$
$`x{\displaystyle \frac{r}{R}}`$
$`\beta {\displaystyle \frac{N}{RT}}`$ (1)
$`uv\sqrt{{\displaystyle \frac{R}{N}}}`$
$`\tau t{\displaystyle \frac{N^{1/2}}{R^{3/2}}}.`$
The unit of time, often referred to as crossing time, $`[t]=\frac{R^{3/2}}{N^{1/2}}`$, is obtained by dividing the unit of length $`R`$ by the unit of velocity $`\sqrt{N/R}`$. This unit of time is also proportional to the period of plasma oscillations in a medium with the charge concentration $`N/R^3`$. As this time unit has purely kinematic origin, we do not expect the evolution of systems having different $`N`$ and $`R`$ to be universal in time $`\tau `$. The evolution, assuming that it is collisional, is expected to be universal in the relaxation time $`\tau _r=\tau \frac{\mathrm{ln}N}{N}`$ bt , where the factor $`N/\mathrm{ln}N`$ is proportional to the number of crossings a typical particle needs to change its velocity by a factor of 2 through weak Coulomb scattering events.

The soft core radius $`x_0=r_0/R=5\times 10^3`$ is chosen to be well below the critical value $`x_{gr}0.021`$, above which the collapse-explosion transition is replaced by a normal first-order phase transition chi .

The MF theory of the system is described in detail in Ref. usg . The equilibrium velocity distribution is Maxwellian and isothermal, while the equilibrium (saddle point) density profile $`\rho (x)`$ corresponding a stable or a metastable state, is a spherically symmetric solution of the integral equation (II). This equation replaces the Poisson-Boltzmann differential equation for the self-consistent potential (see, for example, Ref. pr ) since the interparticle interaction considered here is not pure Coulombic.

$`\rho (𝐱)=\rho _0F[\rho (.),𝐱]`$
$`F[\rho (.),𝐱]=\mathrm{exp}\left[\beta {\displaystyle }{\displaystyle \frac{\rho (𝐱^{})}{\sqrt{(𝐱𝐱^{})^2+x_0^2}}}d^3𝐱^{}\right]`$
$`\beta ={\displaystyle \frac{3}{2}}\left[ϵ+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\rho (𝐱_1)\rho (𝐱_2)}{\sqrt{(𝐱_1𝐱_2)^2+x_0^2}}d^3𝐱_1d^3𝐱_2}\right]^1`$ (2)
$`\rho _0=\left\{{\displaystyle }F[\rho _s(.),𝐱]d^3𝐱\right\}^1`$
The equilibrium density profile $`\rho (x)`$ obtained from this equation is then used to calculate the entropy and the pressure.

The MF phase diagram of the system is presented in Fig. 1. The collapse and explosion energies are $`ϵ_{coll}0.339`$ and $`ϵ_{expl}0.267`$. Examples of the uniform and the core-halo density profiles for $`ϵ=ϵ_{coll}`$ are shown in Fig. 2.

In the MD simulations we consider a system consisting of $`N=125`$ – 500 particles in a spherical container of the radius $`R=1`$. All interparticle forces are calculated directly at each time step $`dt`$. This is done to avoid any mean-field-like effects inevitably present in any truncated multipole or Fourier potential expansion. The particle velocities and coordinates are updated according to the velocity-Verlet algorithm which is symplectic and reversible.

The system is initialized by randomly distributing particles according to a spherically symmetric density profile; typically the appropriate MF density profile $`\rho (x)`$ was used. The potential energy ($`U`$) of the initial configuration is calculated, and the target kinetic energy $`E_k=EU`$ is determined. The particle velocities are randomly generated from some (usually Maxwell) distribution with the appropriate square average. Finally, the deviation of the total energy from its target value, caused by a stochasticity in velocity assignment, is determined, and the velocities are rescaled to fine-tune the total energy. Due to the isotropicity of the random velocity assignment, we have always obtained the states with sufficiently low total angular momentum which collapsed to single-core states rather than to binaries grv .

To implement the reflective boundary condition, at each time step the normal components $`v_{}`$ of the velocities of all particles which had escaped from the container were reversed. Values of the normal components were stored to evaluate the pressure on the wall $`P`$.

$$P(t)=\frac{_{t^{}=tt^{\prime \prime }/2}^{t^{}=t+t^{\prime \prime }/2}v_{}(t^{})}{2\pi R^2t^{\prime \prime }}$$
(3)
During each simulation run we measured such characteristics as the kinetic energy $`ϵ_{kin}=\frac{3}{2\beta }`$, the virial variable $`\sigma `$ (dimension of energy) quantifying deviations from the virial theorem, $`\sigma ϵ+ϵ_{kin}\frac{3PVR}{N^2}`$ (where $`P`$ is the pressure on the wall, $`V=4\pi R^3/3`$ is volume of the container, and the factor $`N^2/R`$ rescales the volume-pressure term to the unit of energy introduced in (II)), ratio of the velocity moments, $`\frac{19v^2^2}{5v^4}`$ (which should be 1 for the Maxwell distribution), and the number of particles in the core $`N_c`$ of the prescribed radius $`x_c`$. For the last measurement we count the number of particles $`N_i`$ that are within $`x_c`$ from the $`i`$th particle and find the particle which has the largest $`N_i`$.

In addition, we measured the histograms of the velocity distribution and radial distribution functions, $`W(u)`$ and $`C(x)`$, respectively. The latter was defined as the number of particles in the spherical layer of the radius $`x`$ around each particle, normalized by the volume of such layer, disregarding the nonuniformity of the system and the boundary effect.

The measurements of the ”scalar” quantities such as energy, kinetic energy, pressure, and velocity distribution moments were taken in time intervals $`\tau _{meas}`$, which were selected sufficiently long to avoid measuring the unchanged configuration repeatedly and sufficiently short not to miss the important details of the system evolution. We usually pick $`\tau _{meas}`$ of the order of the uniform density sphere crossing time $`\tau _{cross}^u=\pi `$, which is a half period of the oscillation of a particle released with zero velocity at the container wall. The histogram data, such as the velocity distribution and the radial distribution functions, was incremented at each $`\tau _{meas}`$ and accumulated over a longer time period $`\tau _{hist}`$, $`\tau _{hist}10`$$`10^3\times \tau _{meas}`$.

Our attempts to resolve the high-density part of the radial density profile of the system turned out to be fruitless due to the strong fluctuations in the positions of this part. This fluctuations result in smearing the central peak in both core-halo and low energy uniform states. Considering the center of mass system of reference does not resolve this difficulty, as, despite being dense, the core typically contains only 10 – 20% of the total system mass (see below) and the positions of the core and the center of mass of the system do not usually coincide.

To control the quality of the simulation, we monitored the total energy $`ϵ`$ and the total angular momentum $`L`$. We selected timestep $`dt`$ small enough to keep the total energy variation within 0.01% of its initial value, usually we used $`dt=10^5`$, or in rescaled units, $`d\tau 10^4`$. For such time steps, the relative deviation of the angular momentum was within $`10^{14}`$.

All the measurements below are presented in the rescaled dimensionless units as defined in Eq. (II).

## III uniform and core-halo equilibrium states: comparison to the MF

To check our simulation procedure and possibly resolve the apparent contradiction between the MF and the particle simulation results pet , we first considered the system in what we expected to be a stable or a metastable states far away from a transition point. Since we were interested in the equilibrium properties, we were initiating the MD systems according the corresponding MF predictions. It meant that the density profiles were seeded according to the MF profiles and the velocities were assigned according to the Maxwell distribution. We observed that the MF density initiation virtually eliminates the transitory period, while the method of velocity assignment was practically unimportant, provided that it gave the correct value for the total kinetic energy. For example, it takes a system initialized with a flat $`W(𝐮)=const`$ about $`\tau \tau _r`$ to evolve to the Maxwell distribution.

A typical plot of the steady state time dependence of the kinetic energy, virial variable, and the total energy is presented in Fig. 3.

The comparison between the MD measurements and the MF results for the uniform and core-halo states is presented in Tables LABEL:tab\_un and LABEL:tab\_ch and reveals a perfect agreement between these two sets of data.

To obtain the expression for the MF virial variable

$$\sigma _{MF}=ϵ+ϵ_{kin}(18\pi \rho (1)/3)$$
(4)
we write for the pressure at the container wall $`P=2\rho (x=1)ϵ_{kin}/3`$ implying that the system is isothermal. Since the interparticle potential is not pure Coulombic, the virial variable is non-zero. The difference is especially prominent for the core-halo states where more particles ”probe” the short-range part of the potential.

To evaluate the core radius and number of core particles of the core-halo system, we considered an integrated MF density profile, $`f(x)=_0^x4\pi y^2\rho (y)𝑑y`$ (see Fig. 4).

As it follows from the Figure, the MF core-halo state indeed contains a distinct core with a sharp boundary of the radius $`x_c10^2`$ relatively insensitive to the energy in the range we considered, $`|ϵ|<0.5`$. Using this MF core radius, we located cores in the MD core-halo systems which contained very similar to the MF cores number of particles (see Table LABEL:tab\_ch). Using smaller core radius resulted in significant reduction in the number of observed core particles. A reasonably small over-estimation of the core radius did not affect the results of the MD measurements: we observed that even in the sphere twice the core radius the number of particles is only marginally (at most by 8%) larger than in the core.

To check if the system has more than one core, we performed search for the second-largest core of the same radius $`x_c`$. We looked for a largest group of particles which are within $`x_c`$ from a single particle with none of these particles belonging to the first, largest core. We never observed the second-largest core containing more than 2 particles; most of the time it contained only a single one.

In Fig. 5 we present the MD velocity distribution functions $`W(u)`$ for a core-halo and uniform states; shown $`W(u)`$ confirm the MF prediction for the Maxwellian form of these distributions.

As we mentioned in the previous section, we were unable to resolve the high-density part of the radial density profile due to the core motion. However, an indirect comparison between the radial distribution of particles in the MF and MD was made using the radial distribution function. The MF radial distribution function $`C_{MF}(x)`$ was computed as

$$C_{MF}(x)=\frac{1}{4\pi x^2}\rho (𝐱^{})\rho (𝐱+𝐱^{})𝑑𝐱^{}.$$
(5)
The good agreement between the MF and the MD radial distribution functions is illustrated in Fig. 6. This indicates that the mutual distribution of particles is correctly predicted by the MF theory.

To summarize, for all the quantities considered, we observed no systematic deviations between the MF theory and the MD data.

## IV Collapse

According to the MF theory, if the energy of the uniform state becomes lower than $`ϵ_{coll}0.339`$, the system should undergo a collapse to a core-halo state. To study the collapse, we considered several uniform systems with the energies ranging between $`ϵ=0.5`$ and $`ϵ=0.3`$. The systems were initialized according to the MF density distributions. For systems with $`ϵ<ϵ_{coll}`$ the particles were distributed according to the MF density profile for $`ϵ_{coll}`$.

In a perfect agreement with the MF theory, a uniform state with $`ϵ<ϵ_{coll}`$ undergoes a gradual transition to a core-halo state with a typical timescale of $`\tau _{coll}10^4`$ for $`N=125`$ – 250 particles. An example of the time dependence of the kinetic energy and the virial variable for a collapsing system is shown in Fig. 7.

We observe that if the number of particles is increased but the rescaled energy $`ϵ`$ is kept fixed, it takes generally longer time for the collapse to be complete. Our results (see Fig. 8) qualitatively confirm that the characteristic time for the full collapse scales as $`\tau _r`$ bt . A quantitative study of the dependence of the collapse dynamics on the number of particles requires much faster simulation code, however.

In the above examples, the energy was set to $`ϵ=0.5`$ which is well below $`ϵ_{coll}0.339`$, and as a consequence the collapse started immediately at $`\tau =0`$ in all simulation runs. If the system energy is $`ϵ_{coll}`$, the noticeable increase in kinetic energy and decrease of the virial variable, characteristic for collapse, start not exactly at $`\tau =0`$ but with a small delay (Fig. 9) which varies from run to run from almost zero to about $`\tau 1500`$. This indicates that the MD system is able to overcome the metastability at or near $`ϵ_{coll}`$. The observed uncertainty is likely due to the relatively small number of particles.

As we increase the energy above $`ϵ_{coll}`$, the stability of the uniform state increases which results in a longer lifetime of such state with respect to collapse. In Fig. 10, an evolution of a system with $`ϵ=0.3`$ is shown. The system stays in the uniform state for about $`\delta \tau 5000`$ before the collapse starts, after which the evolution proceeds qualitatively similar to the collapses in systems with lower energies.

To compare the temporal evolution of the kinetic energy, virial variable, and the number of core particles, the relative variables $`ϵ_{kin}^{}(\tau )`$, $`\sigma _{kin}^{}(\tau )`$, and $`N_{core}^{}(\tau )`$, all defined as $`ϵ_{kin}^{}(t)=[ϵ_{kin}(t)ϵ_{kin}(u)]/[ϵ_{kin}(ch)ϵ_{kin}(u)]`$, are plotted in Fig.11. The values $`ϵ_{kin}(u)`$ and $`ϵ_{kin}(ch)`$ correspond to the uniform and core-halo states in equilibrium.

Figure 11 indicates that during the initial stages of collapse the core grows faster than the kinetic energy and the virial variable. In addition, one can notice large reversible fluctuations in the number of core particles (the core grows up to 12% of its final value and then disappears) that are not matched by comparable scale fluctuations in the kinetic energy or virial variable. All these observations suggest that the density evolution causing the core formation plays the leading role in the process of collapse while the relaxation of kinetic energy follows. Once the collapse has started, the core grows to about a half of its final size in only $`\delta \tau _{core}10^3`$ for systems with $`N=125`$ – 500 particles, while the changes in kinetic energy during this interval of time are small. After this rapid initial stage the system relaxes more slowly, and finally after $`\tau _{coll}10^5`$ reaches the equilibrium core-halo state. Our observations strongly suggest that the growth of the core takes place through a sequential absorption of single particles rather than through hierarchical merging of smaller cores: We never detected other cores containing more than two particles. Although the kinetic energy relaxation trails behind the the core formation, the velocity distribution function remains Maxwellian throughout the whole evolution with the temperature corresponding to the corresponding value of the kinetic energy. This is caused by the fast velocity relaxation ($`\tau _{vel}1`$) as discussed in the previous section.

## V Explosion

It is natural to assume that if a system exhibits a collapse, it should also exhibit an explosion which is the reverse to the collapse transition. According to the MF theory, such explosion should take place when the core-halo state becomes unstable, i.e., when $`ϵϵ_{expl}0.267`$. To check this prediction, we initialized the MD system according to the MF equilibrium core-halo state and followed its evolution. As in the study of the collapse, for initial states with $`ϵ>ϵ_{expl}`$ we used the MF density profiles of the highest energy locally stable state, i.e. of the state with $`ϵ=ϵ_{expl}`$.

We observe that a system with sufficiently high energy, such as $`ϵ=0.5`$ in Fig. 12 or $`ϵ=0.4`$ in Fig. 13, indeed undergoes an explosion which brings it to the uniform equilibrium state. During such an explosion, the state variables such as kinetic energy and virial variable continuously change from their equilibrium core-halo state values to the uniform state ones, and the core gradually sheds particles until only one particle is left.

The main features of an explosion (Figs. 12 and 13) resemble those of a time-reversed collapse. The kinetic energy evolves relatively uniformly, while the number of core particles changes only slightly during the first stages of evolution and rapidly decreases at the final stages. In the example presented in Fig. 12, the explosion is complete after the time $`t_{expl}15000`$, which is noticeably less than the time for a collapse $`t_{coll}10^5`$ (see Fig. 8) for a system having the same number of particles ($`N=250`$). However, the latter is rather vaguely defined due to larger fluctuations in a core-halo than in a uniform state.

Similarly to a collapse, the system remains thermalized in the velocity space during an explosion. The velocity distribution remains Maxwellian throughout the evolution with the temperature corresponding to the current value of kinetic energy. As an illustration, Fig. 14 shows the ratio of the moments of velocity distribution, $`19v^2^2/5v^4`$, which should be 1 for a Gaussian distribution.

However, as it is evident from a comparison between the Figs. 12 and 13, as $`ϵ`$ gets closer to $`ϵ_{expl}`$ the explosion takes longer to initiate. We have observed, that even for $`ϵ=0.3`$, which is noticeably larger than $`ϵ_{expl}0.267`$, the explosion does not happen during the first $`\tau =30000`$ of evolution. This suggests that either the MF value for $`ϵ_{expl}`$ is incorrect, or during the incitation of the system we somehow prepare the system not exactly in the equilibrium (metastable) core-halo state. If the latter is the case, a deviation from the equilibrium most probably takes place in the core, as because of its compactness, its equilibration with the rest of the system may take a rather long time. Using the current MD setup, we were unable to determine a reason for this apparent discrepancy.

## VI CONCLUSION

In the previous sections we have presented the following molecular dynamics results for the self-attracting systems with soft Coulomb potential:

* A collapse from a uniform to a core-halo state was observed. The timescale for the collapse in systems consisting of 125 – 500 particles is of order of $`10^5`$ crossing times and is by the same factor longer than the timescale of the velocity relaxation. The collapse starts with a fast growth of a core via absorption of single particles and continues with more gradual relaxation towards an equilibrium core-halo state. Metastable states exhibit a finite lifetime before collapsing.
* A reverse to collapse, i.e., an explosion transition from a core-halo to a

uniform state was observed. The explosion time is considerably shorter than the collapse time, being of the order of $`10^4`$ crossing times (125 – 500 particles). An explosion resembles a time-reversed collapse; the core decrease, which happens by shedding individual particles, is trailing the kinetic energy evolution till the last stages, when the core rapidly disappears.
* Such molecular dynamics characteristics of the equilibrium or the metastable uniform and core-halo states as kinetic energy, wall pressure, number of core particles, particle-particle radial distribution function and velocity distribution function, are found to be equal within the statistical uncertainty of the molecular dynamics measurements to the corresponding mean field predictions.

The long collapse time observed in our simulations appears to be an explanation for the apparent discrepancy between the phase diagram presented in pet  and the mean field phase diagram (see, for example, chi ). The relaxation time allowed in pet  before the measurements of what was considered to be a steady state, $`t_{rel}=3N/|EN|^{3/2}`$, which is apparently equivalent to $`\tau _{rel}<1`$, is by far insufficient for a system to collapse. Therefore, the discontinuities in caloric curves $`\beta `$ vs $`ϵ`$, typical for collapse and explosion gravitational transitions, were not observed in pet .

Although we considered systems only with the soft Coulomb potential, we speculate that a likewise similarity between the mean field and molecular dynamics equilibrium properties of the core-halo state exists for all ”soft” long-range (like a Fourier-truncated Coulomb) potentials. This is so because all soft potentials are effectively longer-ranged than the bare Coulomb one. However, the core-halo state in the system with a ”harder” short-range cutoff may have completely different properties from the one considered above, and its mean field theory may be inadequate. As for the uniform states, their properties are virtually independent on the nature of the cutoff (see, for example, chi ) and their mean-field description is universally correct.

The main goal in the paper was to check the existence of collapses and explosions and the validity of the mean field data for the self-gravitating systems with short-range cutoff. For this goal one or few molecular dynamics runs for each considered system were sufficient. However, to be able to study the dynamical features of collapses and explosion in more detail and to compare the simulations results to various theoretical models, one needs to study the relaxation averaged over many initial configuration. For example, an interesting question is whether a collapse (or an explosion) indeed consists of two stages; the first fast stage of collisionless ”violent relaxation” with particle number-independent rate, and the slower second stage of soft collisional relaxation with characteristic time $`\tau _r`$ (see, for example, bt  and references therein). Another important question is to resolve the apparent contradiction between the mean field prediction for $`ϵ_{expl}`$ and the molecular dynamics observations, outlined at the end of the previous section. Such studies require a more efficient computation code. The main improvement possibly coming from a better force calculator that may include various mean-field-like potential expansions, which are qualitatively justified by this study. We leave this for the future research.

## VII acknowledgments

The authors are thankful to P.-H. Chavanis and E. G. D. Cohen for helpful and inspiring discussions and gratefully acknowledge the support of Chilean FONDECYT under grants 1020052 and 7020052. M. K. would like to thank the Department of Physics at Universidad de Santiago for warm hospitality.