text:
# Diffusion in a Granular Fluid - Simulation

## I Introduction

Attempts to describe granular media in terms of a more fundamental underlying statistical mechanics have met with considerable success. As a prototype for this approach, in the preceding paper dbyl01  standard linear response methods from normal fluids have been applied, mutatis mutandis, to the case of an impurity particle diffusing in an isolated one component fluid of smooth inelastic hard spheres ($`d=3`$) or disks ($`d=2`$) of diameter $`\sigma `$. The collisions are characterized by a coefficient of normal restitution $`\alpha `$. The isolated system (or with periodic boundary conditions), is not in the typical Gibbs state as for elastic collisions, but rather in a time dependent homogeneous cooling state (HCS). It has been shown in dbyl01 , and will be elaborated further here, that this time dependent HCS can be exactly transformed to a stationary state description. In the present work, this stationary state description is evaluated by molecular dynamics (MD) simulation to measure the mean square displacement of the impurity, its velocity autocorrelation function, and the resulting diffusion coefficient defined in terms of a formal Einstein or Green-Kubo relation, respectively. For practical purposes, attention is restricted to the case of self-diffusion, where the mechanical properties of the impurity are the same as those of the fluid particles.

The velocity autocorrelation function was approximated in dbyl01  using a cumulant expansion and also by means of kinetic theory methods. In the usual first Sonine polynomial expansion for the latter, this leads to a simple exponential decay in an appropriate dimensionless time scale defined below. For the numerical results considered here, one further approximation is made in evaluating the decay rate, namely two-particle velocity correlations are neglected. In this case, the results are the same as those obtained from the Enskog kinetic theory. Earlier studies of self-diffusion BRCyG00  show excellent agreement of this theory with MD results at very low number density $`n`$, over a wide range of values for $`\alpha `$. The present work extends that study to higher densities. More specifically, three-dimensional systems with densities in the interval $`0.1n^{}n\sigma ^30.75`$, for $`0.5\alpha 1`$, will be considered.

The HCS is known to be unstable under long wavelength perturbations or fluctuations GyZ93 . To avoid this problem, the system size in most cases considered is chosen to be smaller than the critical wavelength. The latter is a function of the density and inelasticity, decreasing with increasing values of each. Consequently, at the highest densities and smallest $`\alpha `$ a small system of $`108`$ particles is required. At moderate densities and inelasticities, a system size of $`512`$ particles is used. The agreement of theory and simulation is found to be quite good for all $`\alpha `$ at $`n^{}0.25`$ and for all densities at $`\alpha 0.9`$. At higher densities and smaller $`\alpha `$, significant discrepancies occur, with the diffusion coefficient obtained from MD being almost an order of magnitude greater than the theoretical estimate at $`n^{}=0.75`$ and $`\alpha =0.5`$. The failure of the Enskog theory at high densities is well-known for the case of elastic collisions, due to cage effects and correlated binary collisions. However, in that case the diffusion coefficient from MD is smaller than the theoretical prediction. Therefore, it is clear that a quite different density mechanism is effective for inelastic collisions.

It is reasonable to expect that the underlying instability is responsible in part for the above discrepancies. Although the instability is, in principle, avoided by the use of small systems, the latter is prone to large fluctuations. These are not quantitatively significant at $`\alpha =1`$ but grow rapidly with decreasing $`\alpha `$, as can be seen in the noise level for the kinetic temperature. This will be illustrated later on in Fig. 1. A related issue is the role of the instability in affecting the usual mode coupling corrections to the Enskog theory. Such mode coupling terms arise from correlations of spontaneous fluctuations in the long wavelength hydrodynamic modes, which are significantly modified by the potential instability. Enhanced fluctuations will be addressed in the qualitative analysis of the data below.

The plan of the paper is as follows. In Sec. II, the relevant results in the preceding paper are shortly summarized, and a new scale transformation, useful for practical purposes, is introduced. The equations of motion of the new phase-space variables lead in a direct way to a steady-state simulation method for the HCS. Starting from an arbitrary initial condition, the system rapidly approaches a steady state, whose properties are simply related to those of the HCS. It is important to stress that this relationship is not approximate, but an exact consequence of a change of variables. In particular, the method to measure the mean square displacement of a tagged particle as well as its velocity autocorrelation function is discussed in detail. Moreover, both quantities are shown to lead to equivalent results for the self-diffusion coefficients. Although in practice some discrepancies appear in the numerical results obtained by the two procedures, their origin is well understood.

Comparison of theory and simulation is also started in Sec. II, and completed in Sec. III. As mentioned above, the agreement is fairly good at low densities, but relevant effects, which are not taken into account by the theory developed in dbyl01 , show up even at moderate density. The peculiar nature of these effects for inelastic collisions is discussed. Numerical evidence is provided indicating that the underlying hydrodynamic instability associated to the shear mode plays an essential role in the “anomalous” behavior of the self-diffusion coefficient. Nevertheless, such changes in the diffusion constant do not compromise the existence of the diffusion process, which is confirmed for such conditions by simulation results for the mean square displacement. Finally, Sec. IV contains a short summary of the results and also some indications of the possible extensions of the reported work.

## II Steady-state simulation method

To investigate the nature of diffusion in the HCS and to test the theoretical results presented in the preceding paper dbyl01 , a more fundamental description via MD simulation will be considered. The HCS distribution function for a fluid of $`N+1`$ particles has the scaling form

$$\rho _{hcs}(\mathrm{\Gamma },t)=\left[\mathrm{}v(t)\right]^{d(N+1)}\rho _{hcs}^{}\left(\{𝐪_{ij}/\mathrm{},𝐯_i/v(t)\}\right).$$
(1)
Here and in the following the notation is the same as in ref. dbyl01 . In particular, $`v(t)`$ is the thermal velocity defined in the usual way (with Boltzmann’s constant set equal to one) and $`\mathrm{}`$ is proportional to the mean free path of the gas. The time dependence of the above distribution is due to collisional cooling and is determined from

$$_tv(t)=\frac{1}{2}\zeta (t)v(t),$$
(2)
where $`\zeta (t)`$ is the cooling rate. A direct simulation of the cooling fluid, as described by the Liouville dynamics in the actual phase variables, is difficult, since the rapid cooling of the fluid leads to numerical inaccuracies very soon. One method of dealing with this would be to periodically redefine the time-scale of the simulation, so that the typical particle velocity remains of the order of unity. However, if the rescaling is only used to control the numerical stability of the simulation, the state being simulated would nevertheless be time-dependent. For this reason, some simulation studies employ different type of thermostats, such as an externally imposed Brownian force or thermal boundaries, in order to generate a steady state. While these methods provide more or less realistic models of various experimental procedures, they obviously probe a state which is in some way related to, but not identical with, the homogeneous cooling state.

Soto et al. SyM01 , noting the fact that there is no intrinsic time scale in the dissipative hard sphere model, have proposed to rescale all particle velocities after every collision, thus establishing a steady state similar to that described in dbyl01 . However, this procedure has the effect of replacing the binary collision dynamics of the dissipative hard sphere model by an $`N`$-body dynamics, since the result of a collision is to alter all atomic velocities and not just those of the colliding atoms. It is reasonable to imagine that the binary dynamics is recovered in the infinite system limit, but the connection between the HCS and this dynamics is not clear for the case of the small systems considered in most of the MD simulations. Instead, we follow the procedure used in ref. Lutskopre , according to which an exact mapping of the HCS onto a steady state is exploited as the basis of the simulation method. Following the ideas developed in dbyl01 , all velocities are scaled relative to $`v(t)`$ and the new dimensionless time is given by

$$ds(t)=v(t)dt/\mathrm{}.$$
(3)
This time scale is a measure of the average collision number. The corresponding Liouville equation in these variables supports a stationary HCS solution given by $`\rho _{hcs}^{}\left(\{𝐪_{ij}^{},𝐯_i^{}\}\right)`$, where $`𝐪_i^{}=𝐪_i/\mathrm{}`$ and $`𝐯_i^{}=𝐯_i/v(t)`$. Moreover, the average value of any phase function $`A(\{𝐪,𝐯\})`$ is given by Eq.(41) of reference dbyl01 .

$$A;t=𝑑\mathrm{\Gamma }^{}\rho ^{}(\mathrm{\Gamma }^{})A(\{\mathrm{}𝐪_i^{}\left(s\right),v(t)𝐯_i^{}\left(s\right)\}).$$
(4)
The dynamics in the new phase space is obtained from

$$_s𝐪_i^{}\left(s\right)=𝐯_i^{}\left(s\right),_s𝐯_i^{}\left(s\right)=\frac{1}{2}\zeta ^{}𝐯_i^{}\left(s\right)+L^{}𝐯_i^{}\left(s\right),$$
(5)
where $`L^{}`$ is the dimensionless Liouville operator, and $`\zeta ^{}=\mathrm{}\zeta (t)/v(t)`$ is the dimensionless cooling rate. This is the same as the usual dynamics for hard spheres, except that the term proportional to $`\zeta ^{}/2`$ represents an acceleration between collisions which balances the energy lost during collisions, thus enabling a steady state.

Although it is possible to relate average values of the relevant dynamical functions of the original variables with the average values of the same functions of the scaled variables for a general situation, attention will be restricted in the following to the HCS. For homogeneous systems, $`\zeta ^{}`$ becomes time independent. In principle, the scaled dynamics defined by Eq. (5) can be simulated for the properties of interest, without the complications of continuous cooling in real time. In this formulation, the simulation is qualitatively similar to that for elastic collisions, in the sense that the system rapidly approaches a state for which the time-average of the instantaneous scaled temperature is constant, and subsequent determination of average properties is simpler and numerically more accurate. A technical complication is the need to know the exact value of $`\zeta ^{}`$ a priori, which in general is not possible, since it is determined by the original dynamics. Consequently, it is useful to make a second change of scale,

$$𝐪_i^{}=𝐪_i^{},𝐯_i^{}=\left(w^{}/\zeta ^{}\right)𝐯_i^{},s^{}=\left(\zeta ^{}/w^{}\right)s,A^{}\left(s^{}\right)=A\left(s\right),$$
(6)
where $`w^{}`$ is an arbitrary time independent dimensionless frequency. Equation (5) then becomes

$$_s^{}𝐪_i^{}\left(s^{}\right)=𝐯_i^{}\left(s^{}\right),_s^{}𝐯_i^{}\left(s^{}\right)=\frac{1}{2}w^{}𝐯_i^{}\left(s^{}\right)+L^{}𝐯_i^{}\left(s^{}\right),$$
(7)
with $`L^{}`$ the same $`L^{}`$, but with the replacements $`\{𝐪_i^{},𝐯_i^{}\}\{𝐪_i^{},𝐯_i^{}\}`$. This is the form most convenient for MD simulation, with any reasonable choice for $`w^{}`$. More explicitly, the particle dynamics implied by Eq.(7) and implemented in our event-driven simulations consists of an accelerating streaming between collisions

$`_s^{}𝐪_i^{}`$ $`=`$ $`𝐯_i^{},`$ (8)
$`_s^{}𝐯_i^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}w^{}𝐯_i^{},`$
while the effect of the collision of two particles is to alter their relative velocity according to

$$𝐠_{ij}^{}\stackrel{~}{𝐠}_{ij}^{}=𝐠_{ij}^{}\left(1+\alpha \right)\left(\widehat{𝐪}_{ij}^{}𝐠_{ij}^{}\right)\widehat{𝐪}_{ij}^{},$$
(9)
where $`𝐠_{ij}^{}=𝐯_i^{}𝐯_j^{}`$, $`\widehat{𝐪}_{ij}^{}`$ is the unit vector pointing from atom $`i`$ to atom $`j`$ and the center of mass velocity remains the same.

The relationship of $`\zeta ^{}`$ to $`w^{}`$ is determined by the steady state temperature obtained in the simulation. To see this note that while the instantaneous kinetic energy of the entire system is clearly not constant, a corresponding “temperature” defined as $`T^{}=d^1v_1^2;s^{}^{}`$ does approach a constant, as is shown in Appendix A,

$`T^{}(s^{})`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{w^{}}{\zeta ^{}}}\right)^2\left\{1+\left[{\displaystyle \frac{w^{}}{\zeta ^{}\sqrt{2T^{}(0)}}}1\right]e^{w^{}s^{}/2}\right\}^2`$ (10)
$``$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{w^{}}{\zeta ^{}}}\right)^2.`$
Thus $`\zeta ^{}=w^{}/\sqrt{2T^{}(\mathrm{})}`$ which relates quantities measured during the steady state simulation to those calculated in dbyl01 . In summary, simulation of Eqs.(7) is expected to yield a steady state after a short transient period. Subsequently, ensemble averages of properties can be determined as time averages by making the usual assumption of ergodicity. These properties are directly related to those of the HCS by a simple scale transformation as described above.

This steady-state simulation method removes any limitation on the time for which trajectories may be followed, but there are other limitations due to a long wavelength hydrodynamic instabilityGyZ93  for systems with dimension larger than a critical size $`L_c=2\pi \mathrm{}\sqrt{2\eta ^{}/\zeta ^{}}`$ where $`\eta ^{}=\eta (t)/nm\nu (t)\mathrm{}`$ is the shear viscosity. These instabilities occur when the decay rate of a shear mode fluctuation is less than the cooling rate of the thermal velocity so that their size grows relative to the temperature. In the scaled dynamics, the instabilities show up as ordinary unstable modes that grow exponentially with $`s^{}`$, and are therefore easily identified in the simulations Lutskopre . Since these growing modes represent a spontaneous breaking of the (assumed) spatial homogeneity of the system, the impurity particle no longer undergoes simple diffusion and, therefore, its motion is beyond the scope of this paper. The critical size $`L_c`$ is a function of the density and coefficient of restitution, being smaller for high density and small restitution coefficient. In all the cases reported here, the system size is smaller than that required for the instability to show up. For a system of $`108`$ atoms and $`n^{}0.5`$, this provides no limitation on the accessible densities and coefficients of restitution, otherwise only sufficiently low (high) values of the density (coefficient of restitution) can be studied. A different limitation on the simulation is the “inelastic collapse” McNamaraYoung , which occurs for small values of $`\alpha `$ thereby setting a lower limit on the values accessible to simulation. For the system sizes considered in this paper, the above mentioned limits have been characterized in more detail elsewhere Lutskopre .

The three-dimensional simulations reported in the following begin in all cases with an equilibrated fluid subject to elastic collisions (i.e. $`\alpha =1,w^{}=0`$) with $`T^{}=1/2`$. The coefficient of restitution is then set to the desired value and the scaling parameter $`w^{}`$ is set to $`w^{}=\zeta _E^{}\left(\alpha \right)\sqrt{2T^{}(0)}`$ using the Enskog estimate for the cooling rate $`\zeta _E^{}\left(\alpha \right)`$ as given by Eq. (74) in ref. dbyl01 . This choice ensures that the temperature remains of order $`1`$ for all values of $`\alpha `$, for optimal numerical stability, as well as provides a continuous path to the equilibrium system for which $`\alpha =1`$ and $`w^{}=0`$. The simulation is then continued for $`10^7`$ collisions to allow the system to reach the steady state at $`T^{}=w^2/2\zeta ^2`$ where $`\zeta ^{}`$ is the true cooling rate generated by the simulation. Since the largest system considered consists of $`500`$ atoms, this corresponds to at least $`2\times 10^4`$ collisions per atom, and is more than sufficient to reach the steady state. A final simulation of $`10^7`$ additional collisions is then performed, during which all statistical averages of interest are accumulated. Figure 1 shows the behavior of the instantaneous kinetic energy (or temperature) during several typical simulations. In all cases, the instantaneous temperature fluctuates around a stationary average value, with the size of the fluctuations increasing as the coefficient of restitution decreases. The fact that the average value itself varies with $`\alpha `$ is due to (a) the inadequacy of the simple estimate of $`\zeta ^{}`$ we have used for decreasing values of $`\alpha `$ and (b) to the presence of long-lived shear fluctuations, or vortices, which decay on a time scale of thousands of collisions near the instability. As will be discussed below, the latter give rise to large fluctuations in all components of the kinetic contribution to the pressure tensor as well as raising the apparent temperature.

The system is expected to show a typical diffusive behavior under the scaled dynamics after a short transient time, i.e. for $`s^{}>>1`$. Then the diffusion coefficient can be obtained either from the mean squared displacement (msd) using the Einstein relation, Eq. (70) of ref. dbyl01 , or from the Green-Kubo relation, Eq. (69) of ref.dbyl01 . The former is evaluated, for $`d=3`$, from

$$D^{}(s^{})=\frac{\zeta ^{}}{w^{}}D^{}(s)=\frac{1}{6(N+1)}\frac{}{s^{}}\underset{i=0}{\overset{N+1}{}}\left|𝐪_i^{}(s^{})𝐪_i^{}(0)\right|^2.$$
(11)
Figure 2 shows some typical simulation data exhibiting the expected linearity of the mean-square displacement as a function of $`s^{}`$. The Green-Kubo relation expresses the diffusion coefficient in terms of the velocity autocorrelation function (vacf) as

$$D^{}(s^{})=\frac{1}{3}_0^s^{}𝑑s^{}𝐯_0^{}(s^{})𝐯_0^{}^{}.$$
(12)
The vacf is evaluated from

$$<v_\alpha ^{}(s^{})v_\beta ^{}(0)>=\frac{1}{(N+1)j_{\mathrm{max}}}\underset{i=1}{\overset{N+1}{}}\underset{j=0}{\overset{j_{\mathrm{max}}}{}}v_{i\alpha }^{}(s^{}+j\mathrm{\Delta })v_{i\beta }^{}(j\mathrm{\Delta }),$$
(13)
where $`\alpha ,\beta =x,y,z`$, $`\mathrm{\Delta }`$ is the sampling period, and $`j_{\mathrm{max}}`$ is a function of the amount of data available and the time at which the steady state is established. Notice that the stationarity of the scaled dynamics is explicitly used to evaluate the vacf for a given time separation $`s^{}`$, by averaging over many different samplings during the simulation. In the simulations, $`\mathrm{\Delta }`$ is taken to be $`1/4`$ of the Boltzmann mean free time at $`T^{}=1/2`$, i.e. $`\mathrm{\Delta }^{}=\frac{\sqrt{2}}{4}\left(4n\sigma ^3\sqrt{\pi }\right)^1`$. As discussed above, $`T^{}=1/2`$ is not the exact value of the converged temperature, but in formulating $`\mathrm{\Delta }^{}`$ only a reasonable order of magnitude is required. To evaluate the right side of Eq. (13), we follow the procedure of Futrelle and McGinty Futrelle  (see also the discussion in ref. Allen ), according to which the components of the velocity for each atom sampled at intervals of $`\mathrm{\Delta }^{}`$ are stored. When the number of values stored reaches a threshold, $`N_\mathrm{\Delta }`$, the convolution theorem is used to evaluate Eq. (13) by means of the fast Fourier transform, and so to obtain an estimate of the velocity autocorrelation function. This procedure is then repeated until the end of the simulation and the various estimates, typically between tens and hundreds of samples, are averaged to obtain the final estimate of the correlation function. Figure 3 shows some typical correlation functions obtained in this way, using $`N_\mathrm{\Delta }=1024`$. The functions have been normalized in all cases by their initial values, i.e. the function $`C_{vv}^{}`$, as defined in the preceding paper, has been plotted. To get the diffusion constant, the integral of the velocity correlation function is computed using Simpson’s rule and the variance of the pool of estimates is used as the basis for the calculation of the standard error.

Both methods, Einstein relation and Green-Kubo relation, should yield the same values of the coefficient of self-diffusion. Figures 2 and 3 show that the mean square displacement becomes linear in $`s^{}`$ and the velocity autocorrelation function decays to zero for $`s^{}10`$. For larger $`s^{}`$ both forms for the diffusion coefficient approach a constant, and all further discussion is restricted to this diffusion constant. The consistency between the Green-Kubo and Einstein relations is confirmed in Fig. 4, where the estimates obtained by both procedures are compared for a system of 108 atoms and several values of the density and of the coefficient of restitution. While in qualitative agreement, the difference in values obtained by the two methods is in some cases as great as $`10\%`$. It is expected that the results based on the velocity autocorrelation function are more accurate, because the evaluation at each fixed time $`s^{}`$, is based on an average over all times included in the simulation. In contrast, the mean square displacement is obtained from a single evaluation for each $`s^{}`$. This qualitative difference in accuracy has been confirmed by running multiple simulations for a few systems, and observing variations in the coefficient of self-diffusion as large as $`10\%`$ when determined from the msd, whereas that determined from the vacf showed variations of less than $`1\%`$.

## III Comparison between MD and theory

Figure 5 shows the ratio of the observed diffusion constants to the value predicted by the Enskog-level theory, as a function of $`\alpha `$ for densities $`n^{}=0.1,0.25,0.5`$, and $`0.75`$. For the lowest density, theory and simulation are in excellent agreement even at strong dissipation, as expected from earlier comparisons of simulation results and predictions of the Boltzmann-Enskog equation BRCyG00 . As the density increases, deviations from the Enskog prediction are expected, since they are known to occur even at equilibrium ($`\alpha =1`$). In this respect, our results for $`\alpha =1`$ are quite consistent with previous studies of elastic hard spheresByY80 . It is seen in the figure that deviations increase with decreasing $`\alpha `$. For example, at $`n^{}=0.25`$ the deviations are less than $`10\%`$ for $`\alpha 0.9`$, but increase rapidly to the order of $`25\%`$ for $`\alpha =0.5`$. This behavior can also be observed in the normalized vacf itself which, in the Enskog approximation (here and below we refer to the first Sonine approximation to the Enskog theory), is given by $`\mathrm{exp}(\omega _1^{}s^{})`$, with dbyl01

$$\omega _1^{}=\frac{2(1+\alpha )^2}{3}\chi \left(\pi T^{}\right)^{1/2},$$
(14)
where $`\chi `$ is the pair correlation function for two particles at contact. In Fig. 6, the logarithm of the normalized velocity autocorrelation function, $`C_{vv}^{}(s^{})`$, is plotted as a function of $`\omega _1^{}s^{}`$, for the case of elastic collisions, $`\alpha =1`$, and the several densities we have been considering. Here and below the choice has been made, somewhat arbitrarily but consistently, of truncating the vacf at $`\mathrm{ln}C_{vv}^{}(s^{})4`$ in order to eliminate the tails which are dominated by noise. The data confirm findings of earlier studies ByY80  that the Enskog theory gives a good description at low densities, but increasingly underestimates the diffusion constant, due to the neglect of correlated collisions and cage effects, at higher densities with a maximum deviation around $`n^{}=0.5`$. Above this density, the neglected processes begin to cancel one another and near $`n^{}=0.75`$, the diffusion constant crosses the Enskog prediction, and for higher densities is over-estimated by Enskog theory. Since the Enskog theory is exact at short timesdbyl01 , these effects appear as deviations from the simple exponential form of the vacf at longer times.

Now consider $`\alpha <1`$. A similar plot, Fig. 7, for an inelastic system with $`\alpha =0.7`$, again shows good agreement with the Enskog prediction at the lowest density, but it exhibits much larger deviations than in the elastic case at higher densities, except at short times where the Enskog theory is exact, as already pointed out. Interestingly, the data clearly indicate a crossover to a slower constant decay rate at longer times. Of course, this behavior increases the time integral of the velocity autocorrelation function, and it is responsible for the enhancement of the diffusion coefficient seen in Fig. 5, in contrast to the opposite behavior for high densities at $`\alpha =1`$.

To characterize the deviation of the vacf from the simple exponential form, we start with an exact expression for the vacf based on the Zwanzig-Mori formalism (see Appendix B),

$$_s^{}C_{vv}^{}\left(s^{}\right)+\omega _1^{}C_{vv}^{}\left(s^{}\right)+_0^s^{}𝑑s^{}M(s^{}s^{})C_{vv}^{}(s^{})=0.$$
(15)
where $`M(s)`$ is known as the memory function. If $`M(s)`$ is neglected the exponential decay of the above Enskog approximation is recovered, so the memory function incorporates all of the effects neglected in that approximation. A simple ansatz for this function as an exponential is qualitatively successful in modelling the vacf of fluids with elastic collisions (see ByY80  and references therein), and can be expected to work also for inelastic systems, when formulated in the dimensionless time $`s^{}`$. If we substitute $`M(s^{})=M(0)\mathrm{exp}(\lambda s^{})`$ into Eq.(15), and solve for the vacf with the boundary condition $`C_{vv}^{}(0)=1`$, the resulting model is

$$C^{}(s^{})=\frac{\gamma _{}1}{\gamma _{}\gamma _+}\mathrm{exp}\left(\gamma _+\omega _1s^{}\right)+\frac{\gamma _+1}{\gamma _+\gamma _{}}\mathrm{exp}\left(\gamma _{}\omega _1s^{}\right),$$
(16)
where the constants $`\gamma _+`$ and $`\gamma _{}`$ can be related to $`M(0)`$ and $`\lambda `$. Figure 7 includes the result of fitting these two parameters to the data, being evident that the above model is able to capture the crossover from the Enskog behavior at short times to the slower relaxation for longer times.

While the memory function model provides a framework for describing the results of the simulations, it does not explain them, since any effect not captured in the Enskog theory will give rise to a non-zero contribution to the memory function. In a previous study Lutskopre , similar deviations from the Enskog theory were found for the pressure of the system, and evidence was provided there suggesting that the large discrepancies between theory and simulation for dense, dissipative states may be due to the underlying hydrodynamic instability. Additional support for this possibility comes from Fig. 8 which shows the diffusion constants obtained from the simulation data for systems composed by 500 particles compared to those for systems composed by only 108 particles. For a density of $`n^{}=0.25`$, there is a significant further increase of the diffusion constant in the large system relative to the small system as $`\alpha `$ decreases, with the diffusion constant of the 500 atom system growing to almost twice that of the 108 atom system at $`\alpha =0.5`$. At $`n^{}=0.5`$, the enhancement is even larger, but we know Lutskopre  that the system is in the unstable regime, at least for $`\alpha 0.7`$ and possibly for smaller values, so that most of the enhancement is undoubtedly due to the spontaneous formation of shear flow in the system. This suggests that in the smaller systems, although they are stable, there are present large fluctuations characteristic of the instability, such as spontaneous vortices that form and breakup. In this case, the impurity would find itself in a fluctuating local flow field which could enhance the velocity correlations.

To test whether a local flow field plays a relevant role in the self-diffusion process, two different methods were used to calculate the autocorrelation of the impurity velocity relative to the instantaneous local flow field. In the first method, the instantaneous local flow field was calculated by dividing the simulation cell into $`3^3=27`$ cubic subcells. Then each atom’s velocity relative to the instantaneous average velocity of the fluid in the subcell containing it, was used when calculating both the temperature and the vacf. For the 108 particle system, this procedure leads to a lowering of the measured value of the equilibrium diffusion constant by about $`10\%`$, which we believe to be a finite-size effect due to the removal of a significant fraction of the degrees of freedom of the system. For $`n^{}=0.5`$ it also reduces the steady state temperature from nearly twice the initial temperature (see Fig. 1) to about $`1.5T^{}(0)`$ thus demonstrating the size of the contribution of these fluctuations to the temperature. When the same procedure is applied to a system of 500 atoms, the shift in the equilibrium diffusion constant is negligible. The results obtained by this method are presented in Figs. 9-11. In the graphs, the $`\alpha `$-dependant diffusion constant has been scaled by its measured equilibrium value. The conclusion emerging from the figures is that the $`\alpha `$-dependence of the new diffusion constant is substantially closer to the Enskog form. However, this method may be criticized on the grounds that the statistics of the local velocity field are poor, since the calculations only involve a few atoms in each cell. We, therefore, also consider a second method which is specifically designed to eliminate only the part of the local flow due to the longest wavelength fluctuations in the system. A local flow field $`𝐮`$ is defined by summing over the smallest Fourier components compatible with the size and shape of the simulation cell,

$$u_\alpha \left(𝐫\right)=\underset{\left|𝐤\right|>2\pi /L}{}u_\alpha \left(𝐤\right)e^{i𝐤𝐫},$$
(17)
where $`L^3`$ is the volume of the cubic simulation cell. The Fourier components are determined from the instantaneous velocities of the bulk fluid via

$$u_\alpha \left(𝐤\right)=\frac{2}{N+1}\underset{i}{}v_{i\alpha }e^{i𝐤𝐫_i}.$$
(18)
This method is expected to give better statistics than the cell method discussed above. We then use this flow field to define the relative velocity used in the calculation of the temperature and the vacf. Figures 9-11 also show the resulting diffusion constant when this method of subtracting the instantaneous fluctuations is used. From the figures follows that both methods are mutually consistent over the range $`0.5<\alpha <1`$, and that the agreement between the simulation results and the theoretical prediction from Enskog theory is significantly improved. The effects of subtracting the local flow field are much larger for the inelastic case due to the mechanism responsible for shear instability.

## IV Discussion

In this paper, we have continued the discussion of diffusion in a model granular system begun in ref. dbyl01 . The main motivation has been to show that transformation to the steady-state dynamics allows to carry over many standard methods of nonequilibrium statistical mechanics with relatively minor modifications. Here, the exact correspondence between the usual formulation of the dissipative hard-sphere model of granular fluids and the steady-state dynamics was exploited to formulate a particularly convenient simulation method which eliminates the need for additional complications such as exothermic boundary conditions. We have demonstrated that in the steady-state variables, the homogeneous cooling state exhibits standard diffusive behavior such as the linear increase of the mean-squared displacement with time, and the correspondence between the Einstein and Green-Kubo methods of determining the diffusion constant. In addition, the simplest approximation of Enskog kinetic theory provides excellent agreement at low densities for the whole range of inelasticity $`0.5\alpha 1`$. However, although the formalism and concepts of the statistical mechanics of diffusion have been shown to give an adequate description of diffusion in the HCS, large quantitative deviations from the Enskog kinetic theory predictions have been observed even at relatively moderate densities. At the level of the vacf, it was shown that the deviations are principally due to a crossover from the Enskog time dependence, which is exact for short times, to another, slower decay which nevertheless also appears to be exponential in form. The analytic form of the vacf was shown to be well approximated by modeling the memory function of the diffusive process by a simple exponential as has been used in early studies of memory effects in equilibrium fluids ByY80 .

Physically, we have shown by means of constrained simulations that the deviations from the Enskog model are primarily due to the effect of the longest wavelength velocity modes in the system. For sufficiently large systems, the shear mode is linearly unstable (in the steady-state picture) in the classic sense that the time constant associated with its decay goes to zero as the critical value of $`\alpha `$ is approached from above. A linear stability analysis in the steady-state variablesLutskopre  shows that the shear fluctuations of wavevector $`k`$ decay exponentially with a time constant of $`\eta ^{}k^2\frac{1}{2}\zeta ^{}`$,where $`\eta ^{}`$ is the shear viscosity and the stars indicate quantities expressed in the reduced units of dbyl01 . The dominant $`\alpha `$-dependence of this expression comes from $`\zeta ^{}\left(1\alpha ^2\right)`$ and the size of the system enters through the fact that the smallest non-zero value of the wavevector the system can sample is $`k_{\mathrm{min}}^{}=2\pi /L^{}`$ where $`L^{}`$ is the longest dimension of the simulation cell. For fixed values of $`\alpha <1`$, there will always be a critical value of $`L^{}`$ above which the system is unstable to shear fluctuations. For values of $`\alpha `$ above the critical value, the decay of these fluctuations is nevertheless slowed relative to equilibrium and we therefore reason that, even before the onset of the instability and even in systems too small to exhibit the instability, the slowing down of this mode means that long-lived, long-wavelength fluctuations are present. These fluctuations give rise to the observed enhancement of the diffusion constant. This effect is somewhat analogous to what is observed in turbulent systems since the removal of these modes in the calculation of the diffusion constant removed a large part of the discrepancy from the Enskog prediction. The size dependence of the deviations from the Enskog results supports this conclusion, since larger systems have a higher critical value of $`\alpha `$ because the wavevector can take on smaller values. We conclude that a complete theoretical description of diffusion in HCS will require a model for the memory function taking into account the slow relaxation of the shear modes by means, e.g., of a mode-coupling mechanism.

## V Acknowledgments

The research of JWD was supported by National Science Foundation grant PHY 9722133. JJB acknowledges partial support from the Dirección General de Investigación Científica y Técnica (Spain) through Grant No. PB98-1124. JL acknowledges support from the Université Libre de Bruxelles.

## Appendix A Approach to stationarity

The average value of an observable $`A(\mathrm{\Gamma })`$ for a general homogeneous state $`\rho (\mathrm{\Gamma },t)`$ is given by Eq. (4) above which can be written

$$A;t=𝑑\mathrm{\Gamma }^{}\rho ^{}(\mathrm{\Gamma }^{})A(\{\mathrm{}𝐪_i^{}\left(s\right),v(t)𝐯_i^{}\left(s\right)\})$$
(19)
and

$$_s𝐪_i^{}\left(s\right)=𝐯_i^{}\left(s\right),_s𝐯_i^{}\left(s\right)=\frac{1}{2}\zeta _{hcs}^{}𝐯_i^{}\left(s\right)+L^{}𝐯_i^{}\left(s\right).$$
(20)
The subscript on $`\zeta _{hcs}^{}`$ indicates explicitly that it is the dimensionless cooling rate associated with the HCS

$$\zeta _{hcs}^{}=\frac{\mathrm{}\zeta _{hcs}(t)}{v_{hcs}(t)},$$
(21)
where $`\mathrm{}`$ is the mean free path and $`v_{hcs}(t)=\sqrt{2T_{hcs}(t)/m}`$ is the HCS thermal velocity. The temperature $`T_{hcs}(t)`$ obeys the equation

$$_tT_{hcs}(t)=\zeta _{hcs}(T_{hcs}(t))T_{hcs}(t).$$
(22)
Now make the change of variables

$$𝐪_i^{}=𝐪_i^{},𝐯_i^{}=\left(w^{}/\zeta _{hcs}^{}\right)𝐯_i^{},s^{}=\left(\zeta _{hcs}^{}/w^{}\right)s,A^{}\left(s^{}\right)=A\left(s\right),$$
(23)
to get

$$A;t=𝑑\mathrm{\Gamma }^{}\rho ^{}(\mathrm{\Gamma }^{},s^{})A(\{\mathrm{}𝐪_i^{},v(t)\left(\zeta _{hcs}^{}/w^{}\right)𝐯_i^{}\}),$$
(24)
with the definition $`d\mathrm{\Gamma }^{}\rho ^{}(\mathrm{\Gamma }^{},s^{})d\mathrm{\Gamma }^{}\rho ^{}(\mathrm{\Gamma }^{},s)`$. In particular, the kinetic energy is

$`{\displaystyle \frac{1}{2}}mv^2;t`$ $`=`$ $`v_{hcs}^2(t)\left(\zeta _{hcs}^{}/w^{}\right)^2{\displaystyle 𝑑\mathrm{\Gamma }^{}\rho ^{}(\mathrm{\Gamma }^{},s^{})\frac{1}{2}mv^2}`$ (25)
$``$ $`{\displaystyle \frac{1}{2}}mv_{hcs}^2(t)\left(\zeta _{hcs}^{}/w^{}\right)^2v^2;s^{}^{}.`$
The corresponding temperature is

$$T(t)=\frac{2}{d}\frac{1}{2}mv^2;t=2T_{hcs}(t)\left(\frac{\zeta _{hcs}^{}}{w^{}}\right)^2T^{}(s^{})$$
(26)
where $`T^{}(s^{})d^1v^2;s^{}^{}`$. The time derivative of $`T^{}(s^{})`$ is then found to be

$$\left(_s^{}w^{}\right)T^{}(s^{})=w^{}\frac{\zeta (t)}{\zeta _{hcs}(t)}T^{}(s^{}).$$
(27)
This is still valid for a general homogeneous state. Now assume the existence of a scaling solution, which implies $`\zeta (t)/\zeta _{hcs}(t)=\sqrt{T(t)/T_{hcs}(t)}`$. Then

$$\left(_s^{}w^{}\right)T^{}(s^{})=\sqrt{2}\zeta _{hcs}^{}T^{3/2}(s^{}).$$
(28)
The solution is that given by Eq. (10) of the text.

## Appendix B Memory function model

The dimensionless velocity auto correlation function is defined by

$$C_{vv}^{}(s)\psi \left(s\right)\psi ^{},\psi =\frac{v_{0x}^{}}{\sqrt{v_{0x}^2^{}}}.$$
(29)
Using the notation in ref. dbyl01  this can be written in terms of the generators of the dynamics

$$C_{vv}^{}(s)=𝑑\mathrm{\Gamma }^{}\left(e^^{}s\psi \right)\rho _{hcs}^{}\psi =𝑑\mathrm{\Gamma }^{}\psi e^{\overline{}^{}s}\left(\rho _{hcs}^{}\psi \right).$$
(30)
The detailed forms for the linear operators $`^{}`$ and $`\overline{}^{}`$ are given by Eqs. (36) and (42) of ref. dbyl01  but will not be required here. A projection operator is defined by

$$PX=\rho _{hcs}^{}\psi 𝑑\mathrm{\Gamma }^{}\psi X.$$
(31)
It follows then that

$$PX(s)=\rho _{hcs}^{}\psi C_{vv}^{}(s),$$
(32)
with the choice

$$X(s)=e^{\overline{}^{}s}\left(\rho _{hcs}^{}(\mathrm{\Gamma }^{})\psi \right).$$
(33)
The equation of motion for $`X(s)`$ is

$$\left(_s+\overline{}^{}\right)X\left(s\right)=0$$
(34)
and a closed equation for $`PX(s)`$ is obtained by operating on this equation with $`P`$ and $`Q=1P`$ to get the pair of equations

$`\left(_s+P\overline{}^{}P\right)PX\left(s\right)`$ $`=`$ $`P\overline{}^{}QX\left(s\right),`$
$`\left(_s+Q\overline{}^{}Q\right)QX\left(s\right)`$ $`=`$ $`Q\overline{}^{}PX\left(s\right).`$ (35)
Solving formally for $`QX\left(s\right)`$ in the second equation and substituting into the first gives the desired closed equation for $`PX(s)`$Berne

$$\left(_s+P\overline{}^{}P\right)PX\left(s\right)_0^s𝑑s^{}e^{Q\overline{}^{}Q\left(ss^{}\right)}P\overline{}^{}Q\overline{}^{}PX\left(s^{}\right)=0.$$
(36)
Use of (32) gives the corresponding equation for $`C_{vv}^{}(s)`$

$$_sC_{vv}^{}(s)+\omega _1C_{vv}^{}(s)+_0^s𝑑s^{}M\left(ss^{}\right)C_{vv}^{}(s^{})=0,$$
(37)
with the definitions

$`\omega _1`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Gamma }^{}\psi \overline{}^{}\rho _{hcs}^{}\psi }`$ (38)
$`=`$ $`\left(^{}\psi \right)\psi ^{}`$
and

$$M\left(s\right)=𝑑\mathrm{\Gamma }^{}\psi \overline{}^{}e^{Q\overline{}^{}Qs}Q\overline{}^{}\rho _{hcs}^{}\psi \text{.}$$
(39)
The definition for $`\omega _1`$ is the same as that for the Enskog approximation discussed in dbyl01  and approximated in eq. (14) above. Finally, it is clear that structurally identical expressions would have been obtained had we performed this derivation using the scaled variables introduced in Eq. (6) of the text thus justifying the expression in Eq. (15).


format: text

text:
# A Chandra observation of GRO J1744–28: the bursting pulsar in quiescence

## 1 Introduction

X-ray transients sporadically exhibit very bright outbursts during which their X-ray luminosity can be as high as $`10^{3639}`$ erg s-1. However, most of their time they are in their quiescent state during which they emit X-rays only at a level of $`10^{3034}`$ erg s-1. The mechanisms behind this quiescent X-ray emission are still not understood (see, e.g., Menou et al. 1999; Campana & Stella 2000; Bildsten & Rutledge 2002). The most dominant model for the quiescent properties of non-pulsating neutron star transients, which contain a neutron star with a very low magnetic field ($`<10^{10}`$ Gauss), is the one which assumes that the X-rays below a few keV are due to the cooling of the neutron star after the accretion has stopped (see, e.g., Brown, Bildsten, & Rutledge 1998 and references therein). In the quiescent X-ray spectra of several of those systems an extra power law component above a few keV is present (see, e.g., Asai et al. 1996, 1998; Campana et al. 1998a; Rutledge et al. 2001), but the nature of this component is even more unclear (see, e.g., Campana & Stella 2000 for a discussion).

Detailed studies of the quiescent emission from transient X-ray pulsars (with a neutron star magnetic field strength of $`>10^{11}`$ Gauss) have been inhibited by the lack of detected systems. So far, only two X-ray pulsars have been detected in quiescence, A 0535+26 (Negueruela et al. 2000) and 4U 0155+63 (Campana et al. 2001). Their X-ray luminosities were consistent with the predictions of the cooling neutron star model for the non-pulsating systems, suggesting that also this model might apply to the pulsating ones. However, due to the low statistics of the data a detailed testing of the model could not be performed. To get more insight in the quiescent emission of X-ray pulsars those systems have to be observed with higher sensitivity instruments like Chandra or XMM-Newton. Furthermore, more quiescent X-ray pulsars have to be detected to determine if all of those systems are consistent with the cooling model or that some systems might have different properties which would suggest alternative X-ray production mechanisms (e.g., accretion down to the magnetospheric radius; Stella et al. 1994; Corbet 1996). A good candidate for such studies is the “bursting pulsar” GRO J1744–28.

GRO J1744–28 was discovered in December 1995 with the Burst and Transient Source Experiment (BATSE) aboard the Compton Gamma-Ray Observatory (CGRO; Fishman et al. 1995; Kouveliotou et al. 1996). The source exhibited rapidly repeating, very bright X-ray bursts which are likely due to accretion disk instabilities (e.g., Lewin et al. 1996). The source also exhibited pulsed emission with a pulsation frequency of 2.1 Hz (Finger et al. 1996). Its bursting and pulsating nature lead to the source being called the bursting pulsar. So far, two major outbursts have been detected, one starting in December 1995 and the other one a year later in December 1996. The latter one ended in April 1997 (see, e.g., Woods et al. 1999).

Augusteijn et al. (1997) likely detected GRO J1744–28 using a ROSAT observation performed in March 1996. This observation showed a bright source with a luminosity of $`2\times 10^{37}`$ erg s-1 (for an assumed distance111The distance to the source is unknown but the high column density measured by Nishiuchi et al. (1999) and the proximity of the source on the sky to the center of the Galaxy make a large distance likely. of 8 kpc; 0.1–2.4 keV). Archival ROSAT observations did not detect a source on this position with an upper limit on the luminosity of a few times $`10^{33}`$ erg s-1 (0.1–2.4 keV; Augusteijn et al. 1997). This transient nature of the source makes it very likely that indeed the ROSAT source is GRO J1744–28 despite that no bursts or pulsations were detected (note that the large ROSAT upper limit on the pulsations obtained by Augusteijn et al. 1997 was completely consistent with the strength of the pulsations as measured with BATSE and the Rossi X-ray Timing Explorer \[RXTE\]). Localization of GRO J1744–28 by triangulating the data obtained for this source with Ulysses and BATSE (part of the interplanetary network \[IPN\]) resulted in a position (Hurley et al. 2000) which partly overlapped that of the ROSAT error circle, confirming the identification of the ROSAT source with GRO J1744–28. Cole et al. (1997) and Augusteijn et al. (1997) identified a possible infrared counterpart at the edge of the ROSAT error circle, although its enigmatic properties spurred the suggestion that it might be an instrumental artifact (Augusteijn et al. 1997; but see Cole et al. 1997).

Here we present a Chandra/ACIS-I observation of the region containing GRO J1744–28. We discovered a weak source near the center of the ROSAT error circle, which is likely the quiescent counterpart of the ROSAT transient and therefore likely of GRO J1744–28.

## 2 Observation, analysis, and results

GRO J1744–28 was in the field of view ($``$7.1 arcminutes off-axis) during one of the observations which were obtained as part of the Chandra survey of the Galactic Center region (Wang, Gotthelf, & Lang 2002). This particular observation was performed on 2001 July 18 17:37 \- 20:49 UT with an exposure time of $`10.6`$ ksec. The ACIS-I instrument was used during this observation. To limit the telemetry rate only those photons with energy above 1 keV were transmitted to Earth. The data were analysed using the analysing packet CIAO, version 2.2.1, and the threats listed on the CIAO web pages222Available at http://asc.harvard.edu/ciao/.

The resulting image of the region near GRO J1744–28 is shown in Figure 3. We detected one source in the ROSAT error circle (Augusteijn et al. 1997) and the IPN error ellipse (Hurley et al. 2000) of GRO J1744–28. This strongly suggests that indeed we have detected the quiescent counterpart of GRO J1744–28. The best source position (as obtained with the CIAO tool WAVDETECT) is R.A. = 17h 44m 33.09s and Dec. = –28 44 27.0′′ (for J2000.0). Due to the low number of detected photons the statistical error on the source position is 0.3 arcseconds. However, the satellite pointing error is approximately 0.7 arcseconds (1$`\sigma `$; Aldcroft et al. 2000) and dominates the positional inaccuracy. The proposed infrared counterpart (Cole et al. 1997; Augusteijn et al. 1997) is not consistent with our Chandra position.

The elongated structure of the detected source in Figure 3 is due to the extended point-spread-function for a point source approximately 7 arcminutes off-axis. We detected 52$`\pm `$8 counts (corrected for background) from the source position, resulting in a count rate of 4.9$`\pm `$0.7 $`\times 10^3`$ counts s-1. The source spectrum was extracted using a circle with a radius of 10′′ on the source position333For a $`7^{}`$ off axis source an extraction radius of 10′′ ensures that virtually all source photons are extracted and that no contamination of source photons occurs in the background region.. The background data were obtained by using an annulus on the same position with an inner radius of 10′′ and an outer one of 30′′. The data were rebinned using the FTOOLS routine grppha into bins with a minimum of 5 counts per bin. We employed both the $`\chi ^2`$ and CASH (Cash 1979) statistics to fit the data. Both methods give very similar results and we will only discuss the results obtained with the $`\chi ^2`$ method. The obtained spectrum was fitted using XSPEC version 11.1.0 (Arnaud 1996).

The spectrum (Fig. 3) was of poor quality and every single-component model provided an acceptable fit ($`\chi ^2/dof\begin{array}{c}<\\ \end{array}1`$). The column density could not be constrained and was fixed to the value as determined using ASCA data ($`N_\mathrm{H}5.5\times 10^{22}`$ cm-2; Nishiuchi et al. 1999) during times when the source was in outburst. A pure power-law model resulted in a photon index of 2.4$`\pm `$0.7 (all fit parameters are for a 95% confidence level) and a flux of $`5.3\times 10^{13}`$ erg cm-2 s-1 (0.5–10 keV; unabsorbed; all subsequent fluxes are for this energy range and unabsorbed). A blackbody model resulted in a temperature $`kT`$ of 0.8$`\pm `$0.2 keV, with a radius of the emitting region of only 0.2$`{}_{0.9}{}^{}{}_{}{}^{+1.6}`$ km (assuming a distance of 8 kpc), and a flux of $`2.5\times 10^{13}`$ erg cm-2 s-1. When fitting a neutron star hydrogen atmosphere model (Zavlin, Pavlov, & Shibanov 1996) to the data, the radius of the emitting region could not be constrained and was fixed to 10 km. The obtained temperature at infinity $`kT_{\mathrm{}}`$ was 0.4$`\pm `$0.2 keV and the resulting flux is $`3.0\times 10^{13}`$ erg cm-2 s-1. The uncertainties in which spectral model to use resulted in a possible range for the flux of $`2.55.3\times 10^{13}`$ erg cm-2 s-1. If the source is indeed near the Galactic center at a distance of $`8`$ kpc, then the 0.5–10 keV luminosity would be in the range of $`24\times 10^{33}`$ erg s-1.

When the quiescent spectrum of non-pulsating neutron star X-ray transients are fit with a single component model consisting of either a blackbody or a neutron star atmosphere model, the obtained temperatures are usually 0.2–0.3 keV (using a blackbody model; e.g., Bildsten & Rutledge 2002 and references therein) or $`0.1`$ keV (using a neutron star atmosphere model). Our spectrum of GRO J1744–28 is not consistent with such low temperatures ($`\chi ^2/dof>3`$): both the blackbody and the neutron star atmosphere model fall-off more rapidly at energies above 1 keV than our data. However, several of those non-pulsating systems have displayed a composite quiescent X-ray spectrum, with a soft, most likely thermal component below $`1`$ keV and a hard, power-law type component at higher energies (e.g., Asai et al. 1996, 1998). We have fitted our data using such a composite spectrum consisting of either a blackbody or a neutron star atmosphere model for the soft component and a power-law with photon index of 1 or 2 for the hard component. However, due to the low number of photons detected and the high column density towards the source, no constraints could be set on the temperature of the soft component. The data are fully consistent with the temperatures observed for other systems (e.g., a blackbody temperature of 0.2–0.3 keV) and with an upper limit on the thermal flux of $`3\times 10^{13}`$ erg cm-2 s-1 (0.5–10 keV; for a $`kT`$ of 0.3 keV) or $`9\times 10^{13}`$ erg cm-2 s-1 (for a $`kT`$ of 0.2 keV).

The detected source was the only one present in the ROSAT error circle and the IPN ellipse. Hands et al. (2002) reported on XMM-Newton observations performed on the Galactic plane to study the Galactic X-ray point source population. Using their log $`N`$ – log $`S`$ curve and a 2–6 keV flux of $`12\times 10^{13}`$ erg cm-2 s-1 for GRO J1744–28, we estimate that about 10 to 20 sources per square degree can be detected with the same flux. This gives a probability of $`12\times 10^4`$ that a random field source would fall in the surface area traced out by the intersection of the ROSAT error circle and the IPN error ellipse. The source density will likely be higher in the Galactic center region which would increase this probability. But the calculated probability can be used as a first approximation and it indicates that the detected source is likely GRO J1744–28. However, if the detected Chandra source is not GRO J1744–28, then the flux upper limit on this source would be about $`25\times 10^{14}`$ erg cm-2 s-1 (depending on which spectral model is assumed). Our observation had a time resolution of only $`3.2`$ seconds, which did not allow us to search for the 2.1 Hz pulsations. The low number of photons detected did not allow for a stringent conclusion on the possible variability of the source on longer time scales.

## 3 Discussion

We detected a Chandra source in the ROSAT error circle (Augusteijn et al. 1997) and IPN ellipse (Hurley et al. 2000) of GRO J1744–28, making it likely that we have detected the quiescent X-ray counterpart of GRO J1744–28. The proposed infrared counterpart (Cole et al. 1997; Augusteijn et al. 1997) is not consistent with our Chandra position and the nature of this source is unclear (it might be, as suggested by Augusteijn et al. 1997, an artifact, although Cole et al. 1997 could not confirm this). Our quiescent flux of $`2.55.3\times 10^{13}`$ erg cm-2 s-1 (0.5–10 keV; unabsorbed) is consistent with the non-detection of GRO J1744–28 during the archival ROSAT observation reported by Augusteijn et al. (1997; using their count rate upper limit for GRO J1744–28, an 0.5–10 keV unabsorbed flux upper limit range can be derived using PIMMS444Available at http://heasarc.gsfc.nasa.gov/Tools/w3pimms.html of $`710\times 10^{13}`$ erg cm-2 s-1, depending on the spectral model).

One possible mechanism producing the quiescent emission is residual accretion on to the surface of the neutron star. However, Cui (1997) argued (based on the sudden decrease of the pulsation amplitude at times when the 2–60 keV flux of GRO J1744–28, as measured with RXTE, was below $`2\times 10^9`$ erg cm-2 s-1) that GRO J1744–28 is in the “propeller” regime at relatively low fluxes. In this regime, the magnetic field of the neutron star inhibits accretion on to the neutron star surface. Our Chandra flux for GRO J1744–28 is much lower than the critical flux limit given by Cui (1997) strongly suggesting that the source was in the propeller regime during our observation. Therefore, accretion on to the surface of the neutron star is not likely to be the cause behind the quiescent X-ray emission of GRO J1744–28. Similar conclusions were also reached for the quiescent emission from the transient X-ray pulsars A 0535+26 (Negueruela et al. 2000) and 4U 0155+63 (Campana et al. 2001).

A possible mechanism to produce the observed quiescent X-rays when the source is in the propeller regime might be accretion down to the magnetospheric radius $`r_\mathrm{m}`$ (e.g., Stella et al. 1994; Corbet 1996; Campana et al. 1998b), which is approximately give by $`(GM_{\mathrm{ns}})^{1/7}\mu ^{4/7}\dot{M}^{2/7}`$, in which $`M_{\mathrm{ns}}`$ is the neutron star mass, $`\mu `$ the neutron star magnetic moment, and $`\dot{M}`$ the accretion rate. The luminosity produced will be $`L=GM_{\mathrm{ns}}\dot{M}/r_\mathrm{m}`$. Using our measured luminosity and the magnetic field strength found by Cui (1997; $`2\times 10^{11}`$ Gauss), then $`\dot{M}35\times 10^{15}`$ g s-1 and $`r_\mathrm{m}23\times 10^8`$ cm. This radius is larger than the corotation radius $`r_\mathrm{c}=(GM_{\mathrm{ns}})^{1/3}(P_{\mathrm{spin}}/2\pi )^{2/3}`$ (with $`P_{\mathrm{spin}}`$ the neutron star spin period), which is $`10^8`$ cm for GRO J1744–28, and therefore it fulfills the condition that the source should be in the propeller regime. Although a rough estimate555Large uncertainties are present in the exact luminosity of the source (due to uncertainties in the spectral shape and in the distance to the source), its magnetic field strength, and in certain unknown constants (like the accretion efficiency), which have been assumed to be roughly 1. A full detailed discussion of those uncertainties is beyond the scope of this Letter., it indicates that accretion down to the magnetospheric radius might be able to produce the quiescent X-ray luminosity of GRO J1744–28. However, contributions to the X-ray luminosity might also come from several other mechanisms, like accretion onto the neutron star surface (possibly on localized areas such as the magnetic poles) due to the leakage of matter through the magnetospheric barrier.

Another contribution to the quiescent X-ray emission will be the thermal X-ray emission from the neutron star surface which should give a rock bottom lower limit on the X-ray luminosity. Due to the low statistics of our data and the high column density towards GRO J1744–28, we are not able to accurately probe such a thermal component. However, if the temperature of the thermal component is similar to what has been observed for the quiescent spectra of the non-pulsating neutron star transients, then at least the detected emission in our Chandra observation is mostly due to an another component, which has its origin probably in one of the above discussed mechanisms.

Progress in our understanding of the quiescent X-ray emission can be made for GRO J1744–28 by obtaining longer Chandra or XMM-Newton observations of this source. Such observations will be able to constrain its luminosity and spectral shape much better than we are able to due with the present Chandra data and will determine if the quiescent spectrum of GRO J1744–28 is similar to that of the quiescent non-pulsating systems, or that a fundamental difference is present. Detecting the pulsations in observations with sufficient time resolution will also constrain the mechanisms for the quiescent X-rays in this system. Detecting the pulsations in quiescence would also unambiguously proof that the detected source is GRO J1744–28.

Note added in manuscript: After submission of our paper, we became aware of the paper by Daigne et al. (2002) who reported on a XMM-Newton observation of GRO J1744–28 in quiescence. Their reported position and X-ray flux for the source are consistent with ours.

RW was supported by NASA through Chandra Postdoctoral Fellowship grant number PF9-10010 awarded by CXC, which is operated by SAO for NASA under contract NAS8-39073. WQD is supported by the CXC grant SAO-GO1-2150A. We thank Jon Miller for comments on a previous version of this letter.


text:
# I Introduction

## I Introduction

In recent years polyoxometalates which are practical realization of nanomagnets have become the area of intense research because of their enormous variety of structure and fascinating magnetic properties. A particular aesthetic class is that of the ring-shaped iron(III) compounds denoted as molecular ferric wheels. The decanuclear wheel Fe10 synthesized by Lippard et. al.  may be regarded as a prototype of this class. Meanwhile synthesis of ferric wheels with 6, 8, 12 and 18 have also been reported . These materials have dominant anti-ferromagnetic coupling between $`\mathrm{Fe}(\mathrm{III})`$ spins and a singlet ground state. The magnetic properties of such nanoscopic molecules result from the interplay of superexchange interactions between the atomic spins, dipolar coupling of the local moments and on-site anisotropies arising from ligand configurations. The emergence of ferric wheels has led to a renewed interest in the properties of the Heisenberg chain, especially for large spin values. In 1964 Fisher used a classical treatment , in which spin quantization is absent, to study the properties of Heisenberg Hamiltonian and found analytical solutions for the thermodynamic properties of the system. Numerous quantum mechanical calculations have been made for Heisenberg Hamiltonian, mostly for spin-1/2 chains. Calculations for larger $`\mathrm{S}`$ became more interesting after Haldane’s conjecture  regarding the presence of an energy gap in the excitation spectrum from the ground state for integer $`\mathrm{S}`$ chains.

To fit the experimental temperature dependence of the magnetic susceptibility data for $`\mathrm{Fe}_{10}`$ Lippard et. al. adopted a classical spin treatment to obtain the value of exchange interaction parameter $`\mathrm{J}`$. But below 50 K this treatment fails. Ferric wheels with S=5/2, and system size of up to 8 sites have been treated exactly using irreducible tensor operator technique  with the aid of point group symmetry as well as by using the invariance of spin Hamiltonian with regard to interchange of spin sites . Recently there have been theoretical studies to explain the low temperature magnetic susceptibility data of $`\mathrm{Fe}_{10}`$ . The magnetization of ferric wheels exhibit a step like field dependence at low temperature due to the occurrence of field induced ground state level crossing, a spectacular manifestation of quantum size effect in these nanomagnets . There have been NMR and specific heat studies of these ferric wheels to investigate the energy level structure . In appropriate parameter regime these ring systems are considered to be the candidates for the observation of macroscopic quantum phenomena, in the form of quantum coherent tunneling of the Neel vector . To understand these low temperature properties of ferric wheels in detail we need to know the low lying eigenvalue spectrum for these systems. It is also of interest to compare and contrast the zero temperature density of state of S=5/2 rings with that of the exactly known S=1/2 chain. In this paper we have used spin parity together with rotational symmetry of the ferric wheels to obtain the low-lying eigenvalue spectrum for rings with 6, 8, 10 and 12 spin-5/2. The dispersion spectrum reveals interesting features. There have been recent reports of solving the low-lying eigen-spectrum of $`\mathrm{Fe}_{10}`$ using the Density Matrix Renormalization Group . To the best of our knowledge ours is the first study of of ferric wheels using exact methods up to a ring size of 12, S=5/2 iron(III) ions.

We have also studied the quantum spin dynamics of ferric wheels in presence of an alternating field after setting up the Hamiltonian matrix in the subspace of low-lying states. This Hamiltonian includes multipolar spin-spin interaction terms besides a time varying magnetic field (Zeeman term). We have then evolved an initial state, which is taken to be the ground state with a specific value of $`M_S`$ (the $`z`$-component of the total spin) in the absence of the magnetic field, by using the time-dependent formulation of the problem in the restricted subspace. We observe a nontrivial oscillation of the magnetization whose frequency depends on the amplitude of the alternating field. This phenomena is very similar to what has been observed in case of uniaxial magnets .

## II Model Hamiltonian and Computational Details

### A Symmetry Adaptation of Correlated States

The full symmetry of electronic states is of central interest in quantum theory. Since nonrelativistic Hamiltonians $`\widehat{\mathrm{H}}`$ does not depend explicitly on spin, molecular eigenstates can be labeled by the total spin $`\mathrm{S}`$ and the appropriate irreducible representation of the point group. The model Hamiltonian employed in this study is the isotropic exchange Hamiltonian involving exchange interaction between nearest neighbors,

$$\widehat{\mathrm{H}}=\underset{<\mathrm{ij}>}{\mathrm{\Sigma }}\mathrm{J}_{\mathrm{ij}}\widehat{\mathrm{s}}_\mathrm{i}\widehat{\mathrm{s}}_\mathrm{j}$$
(1)
where the exchange interaction $`\mathrm{J}_{\mathrm{ij}}`$ takes the values dictated by experimental studies of structure and magnetic properties. For $`\mathrm{Fe}_6`$ the $`J`$ value depends on the central alkali metal atom : for Na:$`\mathrm{Fe}_6`$ $`J`$=32.77 K whereas for Li:$`\mathrm{Fe}_6`$ $`J`$=20.83 K. In $`\mathrm{Fe}_8`$, $`\mathrm{Fe}_{10}`$ and $`\mathrm{Fe}_{12}`$ $`J`$ is 22 K,15.56 K and 31.97 K respectively .

The total dimensionality of the Fock space of the ferric wheel is given by

$$\mathrm{D}_\mathrm{F}=\underset{1}{\overset{n}{\mathrm{\Pi }}}(2\mathrm{S}_\mathrm{i}+1)$$
(2)
where ’n’ is the total number of spins in the wheel and $`\mathrm{S}_\mathrm{i}`$ is the spin on each ion. In case of $`\mathrm{Fe}_{10}`$, there are 10 iron(III) ions with the dimensionality of the Fock space 60,466,176. In case of $`\mathrm{Fe}_{12}`$, which we have solved exactly this dimension is 2,176,782,336.

Specializing to given total $`M_S`$ leads to Hilbert space dimensionalities which are lower than the Fock space dimensionality. In the case of the $`\mathrm{Fe}_{12}`$ cluster the $`M_S=0`$ space has a dimensionality of about 145 million (144,840,476). The major challenge in exact computation of the eigenvalues, and properties of these spin clusters lies in handling such large bases and the associated matrices. While the dimensions look overwhelming, the matrices that represent the operators in these spaces are rather sparse. Usually, the number of nonzero elements in a row is of the order of the number of exchange constants in the Hamiltonian. This sparseness of the matrices allows one to handle fairly large systems. However, in the case of spin problems, generating the bases states and using the symmetries of the problem is nontrivial.

The isotropic exchange Hamiltonians conserve total spin, S, besides the z-component of the total spin, $`\mathrm{M}_\mathrm{S}`$. Besides these symmetries, the geometry of the cluster also leads to spatial symmetries which can often be exploited. The simplest way of generating bases functions which conserve total spin is the VB method that employs the Rumer-Pauling rule . It is quite easy to generalize the Rumer-Pauling rules to a cluster consisting of objects with different spins to obtain states with desired total spin, S. However, setting-up the Hamiltonian matrix in such a basis can be computationally intensive since the exchange operators operating on a ”legal” VB diagram (diagram that obeys Rumer-Pauling rules) could lead to ”illegal” VB diagrams and resolving these ”illegal” VB diagrams into ”legal” diagrams would present the major bottle-neck. Indeed, the same difficulty is encountered when spatial symmetry operators operate on a VB function. Thus, the extended VB methods are not favored whenever one wishes to apply it to a motley collection of spins or when one wishes to exploit some general spatial symmetries that may exist in the cluster .

It is advantageous to partition the spaces into different total spin spaces because of the usually small energy gaps between total spin states which differ in S by unity. To avoid the difficulties involved in working with total spin eigenfunctions, we exploit parity symmetry in the systems. The parity operation involves changing the z-component of all the spins in the cluster from $`\mathrm{M}_{\mathrm{S}_\mathrm{i}}`$ to $`\mathrm{M}_{\mathrm{S}_\mathrm{i}}`$. There is an associated phase factor with this operation given by $`(1)^{\mathrm{S}_{\mathrm{tot}}+_\mathrm{i}\mathrm{S}_\mathrm{i}}`$. The isotropic exchange operator remains invariant under this operation. If this symmetry is employed in the $`\mathrm{M}_\mathrm{S}=0`$ subspace, the subspace is divided into ”even” and ”odd” parity spaces depending upon the sign of the character under the irreducible representation of the parity group. The space which corresponds to even (odd) total spin we call the even (odd) parity space. Thus, employing parity allows partial spin symmetry adaptation which separates successive total spin spaces, without introducing the complications encountered in the VB bases. However, the VB method can lead to complete factorization of the spin space leading to smaller complete subspaces.

In the ferric wheels, besides spin symmetries, there also exists spatial symmetries. The topology of the exchange interaction leads to a $`\mathrm{C}_\mathrm{n}`$ point group symmetry, where ’n’ is the number of iron ions in the ring. Hence, $`\mathrm{Fe}_6,\mathrm{Fe}_8,\mathrm{Fe}_{10},\mathrm{Fe}_{12}`$ will have C6, C8, C10 and C12 symmetry respectively. It should be mentioned that among these ferric wheels only $`\mathrm{Fe}_{10}`$ is strictly having a ten fold rotational symmetry, rest of them approximately have the above mentioned symmetry. For computational advantage we have assumed the rotational symmetry for all the ferric wheels. This point group appears at first site to present difficulties because the characters in the irreducible representation are in some cases complex which could lead to complex bases functions. This, however, can be avoided by recognizing that in the $`\mathrm{C}_\mathrm{n}`$ group, states with wavevectors ’k’ and ’-k’ are degenerate. We can therefore construct a linear combination of the ’k’ and ’-k’ states which is real. The symmetry representations in the $`\mathrm{C}_\mathrm{n}`$ group would then correspond to the labels $`A`$, $`B`$ and $`E`$, with the characters in the $`E`$ representation given by $`2\mathrm{cos}(rk)`$ under the symmetry operation $`\mathrm{C}_\mathrm{n}^\mathrm{r}`$, with $`\mathrm{k}=\pi /\mathrm{n}`$. The parity operation commutes with the spatial symmetry operations and the full point group of the system would then correspond to the direct product of the two groups. Since both parity and spatial symmetries can be easily incorporated in a constant $`\mathrm{M}_\mathrm{S}`$ basis, we do not encounter the difficulties endemic to the VB theory.

The generation of the complete basis in a given Hilbert space requires a simple representation of a state on the computer. This is achieved by associating with every state a unique integer. In this integer, we associate $`\mathrm{n}_\mathrm{i}`$ bits with spin $`\mathrm{s}_\mathrm{i}`$, such that $`\mathrm{n}_\mathrm{i}`$ is the smallest integer for which $`2^{\mathrm{n}_\mathrm{i}}(2\mathrm{s}_\mathrm{i}+1)`$. In the integer that represents the state of the cluster, we ensure that these $`\mathrm{n}_\mathrm{i}`$ bits do not take values which lead to the $`\mathrm{n}_\mathrm{i}`$ bit integer value exceeding $`(2\mathrm{s}_\mathrm{i}+1)`$. For each of the allowed bit states of the $`\mathrm{n}_\mathrm{i}`$ bit integer, we associate an $`\mathrm{M}_{\mathrm{S}_\mathrm{i}}`$ value between $`\mathrm{s}_\mathrm{i}\mathrm{and}\mathrm{s}_\mathrm{i}`$. For a spin cluster of ’n’ spins, we scan all integers of bit length $`\mathrm{N}=\underset{\mathrm{i}=1}{\overset{n}{\mathrm{\Sigma }}}\mathrm{n}_\mathrm{i}`$ and verify if it represents a basis state with the desired $`\mathrm{M}_\mathrm{S}`$ value. In Fig. 1, we show a few basis functions with specified $`\mathrm{M}_\mathrm{s}`$ value for some typical ferric wheels along with their bit representations and the corresponding integers. Generation of the bases states is usually a very fast step, computationally. Generating the basis as an ordered sequence of integers that represent them also allows for rapid generation of the Hamiltonian matrix elements as will be seen later.

Symmetrization of the basis by incorporating parity and spatial symmetries involves operation on the constant $`\mathrm{M}_\mathrm{S}=0`$ basis by the symmetry operators. Since spatial symmetry operators permute the positions of equivalent spins, every spatial symmetry operator operating on a basis function generates another basis function. Every symmetry operator can be represented by a correspondence vector whose $`\mathrm{i}^{\mathrm{th}}`$ entry gives the state that results from operating on the $`\mathrm{i}^{\mathrm{th}}`$ state by the chosen operator. This is also true for the parity operator, in the $`\mathrm{M}_\mathrm{S}=0`$ subspace.

The first step in constructing symmetry adapted linear combinations is to represent the symmetry operators in the chosen basis as matrices. In our case, the symmetry operators are such that symmetry operation by any operator on a basis state leads to a resultant which is a single basis state. Thus all our symmetry operators can be represented as vectors; the entry in position ’i’ gives the index of the basis function generated by symmetry operation on the basis state ’i’. Since the basis is very large, it is prohibitive to store and manipulate the full basis together with all the associated symmetry vectors. To avoid these difficulties, we have constructed the symmetry matrices in small invariant subspaces. These invariant subspaces are obtained by sequentially choosing a state and operating on it by all the symmetry operators. This gives rise to a set of states on which we again operate by all the symmetry operators, and continue this process until no new basis states are generated. The collection of all these basis states resulting at the end of this process is the invariant subspace. We can set-up symmetry combinations of the basis states in each of the invariant subspaces independently. The symmetry combinations can now be obtained operating on each state by the group theoretic projection operator,

$$\widehat{P}_{\mathrm{\Gamma }_i}=\frac{1}{h}\underset{𝑅}{\mathrm{\Sigma }}\chi _{\mathrm{\Gamma }_i}(R)\widehat{R}$$
(3)
on each of the basis states of the invariant subspace. Here $`\mathrm{\Gamma }_i`$ is the $`i^{th}`$ irreducible representation, $`\widehat{R}`$ is the symmetry operation of the group and $`\chi _{\mathrm{\Gamma }_i}(R)`$ is the character under $`\widehat{R}`$ in the irreducible representation $`\mathrm{\Gamma }_i`$. This process is repeated with the next basis state that has not appeared in any of the invariant subspaces already constructed. The process comes to an end when all the basis states have appeared in any one of the invariant subspaces.

The resulting symmetrized basis is usually over complete in each of the invariant subspaces. The linear dependencies can be eliminated by a Gram-Schmidt orthonormalization procedure. However, in most cases, ensuring that a given basis function does not appear more than once in a symmetrized basis is sufficient to guarantee linear independence and weed out the linearly dependent states. A good check on the procedure is to ensure that the dimensionality of the symmetrized space in the invariant subspace agrees with that calculated from the traces of the reducible representation obtained from the matrices corresponding to the symmetry operators in the chosen invariant subspace. Besides, the sum of the dimensionalities of the symmetrized spaces should correspond to the dimensionality of the unsymmetrized invariant subspace in each of these subspaces.

Generation of the Hamiltonian matrix is rather straight forward and involves operation of the Hamiltonian operator on the symmetry adapted basis. This results in the matrix $`𝐒H`$, where $`𝐒`$ is the symmetrization matrix representing the operator $`\widehat{P}_{\mathrm{\Gamma }_i}`$ and $`𝐇`$ is the matrix whose elements $`h_{ij}`$ are defined by

$$\widehat{\mathrm{H}}|\mathrm{i}>=\underset{j}{\mathrm{\Sigma }}\mathrm{h}_{\mathrm{ij}}|\mathrm{j}>.$$
(4)
The states {i} correspond to the unsymmetrized basis functions. The Hamiltonian matrix in the symmetrized basis is obtained by right multiplying the matrix $`𝐒H`$ by $`𝐒^{}`$. The resulting symmetric Hamiltonian matrix is stored in the sparse matrix form and the matrix eigenvalue problem is solved using the Davidson algorithm .

Computation of the properties is easily done by transforming the eigenstate in the symmetrized basis into that in the unsymmetrized basis. Since the operation by any combination of spin operators on the unsymmetrized basis can be carried out, all relevant static properties in different eigenstates can be obtained in a straightforward manner.

### B Quantum Spin Dynamics

We have studied the dynamics in ferric wheels by setting up the Hamiltonian matrix in the desired $`\mathrm{M}_\mathrm{S}`$ space, which in all cases is restricted to $`\mathrm{M}_\mathrm{S}=0,1`$ and $`2`$. In each subspace we have obtained a few low-lying states using the Davidson algorithm . We have also calculated the spin-spin correlation functions in each of the states. Using the spin-spin correlation functions, we have computed the expectation value of $`\mathrm{S}_{\mathrm{total}}^2`$ operator, from which we have identified the total spin of the state. We observe from the eigenvalue spectrum of all the ferric wheels that ground state and the first excited state are spin singlet (S=0) and triplet (S=1) state respectively, they belong to two different spatial symmetry subspaces. So they will not mix unless there is a perturbation which spoils the $`\mathrm{C}_{10}`$ symmetry of the molecule. We also notice that the first (triplet) and the second (quintet) excited states again fall into different symmetry subspaces. This is true in all the ferric wheels we have studied.

To study quantum dynamics we have considered the following Hamiltonian ,

$$\widehat{\mathrm{H}}=\mathrm{E}_\mathrm{s}\mathrm{D}\widehat{\mathrm{S}}_{\mathrm{z},\mathrm{total}}^2+\mathrm{c}(\widehat{\mathrm{S}}_{\mathrm{x},\mathrm{total}}^2\widehat{\mathrm{S}}_{\mathrm{y},\mathrm{total}}^2)\mathrm{g}\mathrm{h}(\mathrm{t})\widehat{\mathrm{S}}_{\mathrm{z},\mathrm{total}}$$
(5)
Here D is the quadratic anisotropy factor, $`g`$ is the Landé $`g`$-factors for the iron(III) spin respectively, and h(t) is the time-dependent magnetic field, expressed as $`\mathrm{h}(\mathrm{t})=\mathrm{H}_0\mathrm{cos}(\omega \mathrm{t})`$, where $`\omega `$ is the frequency at which field is ramped and $`\mathrm{H}_0`$ is the amplitude of the field. We have chosen D = $`8.8\times 10^3`$ and c = $`10^3`$ (in units of J) in accordance with the experimental values for $`\mathrm{Fe}_{10}`$. We take $`g=2.0`$. The constants $`\mathrm{E}_\mathrm{s}`$ in Eq. (5) correspond to the lowest energies obtained from Eq. (1). The second order anisotropy term allows transitions between states with $`\mathrm{\Delta }\mathrm{M}_\mathrm{S}=\pm 2`$. Both second and third terms in Eq. (5) arises due to magnetocrystalline anisotropy. The exact form of the anisotropy in ferric wheel is not very well established. We have included anisotropy terms only up to second order in the spin variables. The anisotropy in the plane can be formed artificially, e.g. by means of external electric or magnetic fields, pressure or using anisotropic substrate .

The Hamiltonian in Eq. (5) does not allow the mixing of the (singlet) ground state, (triplet) first excited state and (quintet) second excited state in $`\mathrm{Fe}_{10}`$ because of symmetry reason already mentioned earlier. To observe spin dynamics we apply a constant magnetic field that leads to triplet being the ground state of the system. This constant shifting field plays here the role of a chemical potential and is used to provide the optimal conditions for observing the considered quantum spin dynamics. We can represent the external field as

$$\mathrm{B}=\mathrm{B}_0+\mathrm{H}_0\mathrm{cos}(\omega \mathrm{t})$$
(6)
To study the evolution of magnetization as a function of applied oscillatory field we start from a fixed axial field amplitude (in units of $`\mathrm{J}/\mathrm{}`$) and time is increased in small time step of $`\mathrm{\Delta }`$t =0.01 ($`\mathrm{}/\mathrm{J}`$), frequency of the field is kept constant.

We have studied the time evolution of the system by solving the time-dependent Schrödinger equation,

$$\mathrm{i}\mathrm{}\frac{\mathrm{d}\psi }{\mathrm{dt}}=\widehat{\mathrm{H}}(\mathrm{t})\psi .$$
(7)
We assume the system to be in the state with $`S=1,M_S=1`$ at time $`t=0`$. This is the initial state which is time evolved according to the equation

$$\psi (t+\mathrm{\Delta }t)=e^{i\widehat{H}(t+\frac{\mathrm{\Delta }t}{2})\mathrm{\Delta }t/\mathrm{}}\psi (t).$$
(8)
The evolution is carried out by explicit diagonalization of the Hamiltonian matrix $`𝐇(t+\frac{\mathrm{\Delta }t}{2})`$, and using the resulting eigenvalues and eigenvectors to evaluate the matrix of the time evolution operator $`e^{i\widehat{H}(t+\frac{\mathrm{\Delta }t}{2})\mathrm{\Delta }t/\mathrm{}}`$. We set up the Hamiltonian matrix for time evolution in the truncated basis of 3 states corresponding to total spin $`S`$=1. We repeatedly carry out the the time evolution in small time steps of size $`\mathrm{\Delta }t`$ to obtain time evolution over longer periods.

## III Results and Discussion

### A Analysis of Low-lying Spectrum

We have solved the exchange Hamiltonian (Eq. (1)) exactly using the method mentioned earlier to obtain the low-lying eigenvalue spectrum for 6, 8, 10 and 12 site iron(III) rings. We find that ground state, first, second and third excited states are respectively singlet, triplet, quintet and heptet for all ferric wheels. We notice that there is no accidental degeneracy between the energy levels belonging to different symmetry subspaces. Ground state switches between A+ (k =0) and B+ (k = $`\pi `$) subspace for even and odd value of ’m’ respectively, in N=2m ring systems. This was also observed by Mattheiss in case of a spin-1/2 chain and can be understood from Marshall’s sign rule. The gap between the ground state and the first excited state is shown in Fig. 2 as a function of inverse ring size. According to Haldane conjecture, the gap should extrapolate to zero. The extrapolated value, while small is still finite, suggesting that in these rings the finite size effects are still at play at the ring sizes we have studied.

Using the exchange constants estimated for the different ring systems, we estimate the gap between the ground state and the first excited state to be 22.67 K, 11.81 K, and 6.88 K for $`\mathrm{Na}:\mathrm{Fe}_6`$, $`\mathrm{Fe}_8`$, $`\mathrm{Fe}_{10}`$ respectively. Our calculated values compare very well with the experimental values, which are 22.0 K, 12.1 K and 6.45 K for $`\mathrm{Na}:\mathrm{Fe}_6`$, $`\mathrm{Fe}_8`$ and $`\mathrm{Fe}_{10}`$ respectively . This agreement shows that for all practical purposes ferric wheels can safely be assumed as rings neglecting the slight deviation from the exact circular geometry. Our calculated gap for and $`\mathrm{Fe}_{12}`$ is 12.09 K corresponding to the exchange constants predicted from experiments. However, experimental estimate of this gap is lacking.

If we define $`\delta _i`$ to be the energy difference between the $`i`$-th excited state and the ground state, we find from Fig. 3 that the following relationship is satisfied for the ferric wheels,

$$\delta _i=\frac{S_i(S_i+1)}{2}\times E_1$$
(9)
where E1 is the energy gap between the ground state and the first excited state. This indicates that the lowest spin state obeys the Lande interval rule, in agreement with Taft’s conjecture. If we assume $`E_1`$ to be the inverse of the moment of inertia, then the above expression gives the rotational energy of a rigid rotor in a state with quantum number $`S_i`$. Thus, the spin states of ferric wheels can be viewed as the quantized states of a rigid rotor.

In Fig. 4 we have shown the dispersion spectrum of ferric wheels. The value of $`k`$ corresponding to a wavefunction $`\psi `$ can be defined as,

$$T\psi =e^{ik}\psi $$
(10)
where $`T`$ is the translation operator which in case of a ring rotates the ring by one lattice spacing. We have used the spatial $`C_n`$ symmetry of the ferric wheels, which enables us to identify the value of $`k`$ easily for a specified eigenfunction. We find that, in the momentum ($`k`$) sector which contains the ground state, the lowest excitation is to a quintet state. For all the other $`k`$-sectors the lowest excitation is to a triplet state. Previous studies  on antiferromagnetic spin-1/2 Heisenberg chains show that the excitation spectrum is given by $`\mathrm{}\omega =(\pi /2)J|\mathrm{sin}k|`$, where $`k`$ is the wave vector of the excited states measured with respect to that of the ground state. The simple antiferromagnetic spin-wave theory based on the use of the Holstein-Primakoff transformation for each sublattice, leads to excitation spectrum ($`S`$ is the magnitude of the individual spin $`S_i`$),

$$\mathrm{}\omega =2JS|\mathrm{sin}k|.$$
(11)
This relation is supposed to be correct for $`S\mathrm{}`$. We notice that the excitation spectrum for ferric wheels can be fitted to a $`|\mathrm{sin}k|`$ kind of function. Data points in $`k`$=0 or $`k=\pi `$ diverges from the above sinusoidal function. This may be a finite size effect. In the thermodynamic limit of infinite chain length excitation spectrum of ferric wheel will be similar to Eq. (11).

We have also calculated the spin-spin correlation function ($`_{ij}<S_i^zS_j^z>`$) of ferric wheels and Fourier transformed it to find out the structure factor ($`S(q)`$). In Fig 5.(a - d) we have plotted the structure factor as a function of wave vector $`q`$ for different symmetry subspaces. In all the cases $`S(q)`$ shows a peak at $`q=\pi `$. The ground state is a Neél ordered state. It also signifies that ground state will undergo a spontaneous distortion with wave vector $`\pi `$. We find that the weight of $`S(q)`$ is more in case of ground state than excited states.

### B Evolution of Magnetization in presence of an ac field

We have studied the evolution of magnetization in presence of an axial ac magnetic field varying the amplitude of the field. We have kept the frequency of the field fixed at $`\omega `$=10-3. We calculate the magnetization at each time step. When we draw a smooth curve for the time evolution over long time periods, we find a sinusoidal motion,

$$M(t)\mathrm{cos}(\mathrm{\Omega }t)$$
(12)
which can be seen in Fig. 6. Unexpectedly the frequency $`\mathrm{\Omega }`$ of this sinusoidal motion does not correspond to an eigenfrequency of the system or to the period of the external field. When we change the amplitude H0 of the field, the period of the magnetization changes, which is shown in Fig 6. We find a regular dependence of $`\mathrm{\Omega }`$ as a function of H0, as in Fig 7. Frequency of oscillation $`\mathrm{\Omega }`$ becomes very small for some values of H0, again it increases but to a lower value than the previous $`\mathrm{\Omega }_{max}`$. $`\mathrm{\Omega }_{max}`$ values can be fitted to a exponentially decaying curve. This nontrivial resonance can be understood from the viewpoint of Floquet’s theorem . Nonadiabatic transitions are possible whenever system becomes metastable. This resonance can be observed if two states with different magnetizations are nearly degenerate leading to large tunneling cross sections between the states resulting in oscillation of the magnetization. This magnetization oscillation can be related to that of macroscopic quantum coherence, which is predicted by Zvezdin et. al. .

### C Torque Magnetometry

Cornia et. al. have used a novel cantilever torque magnetometry technique to study the spin-state crossover in ferric wheels. The torque T experienced by a magnetically anisotropic substance in a uniform magnetic field B is given by

$$𝐓=𝐌\times 𝐁$$
(13)
where M is the magnetization of the sample. T vanishes when the magnetic field is applied along one of the principle directions $`\widehat{x}`$, $`\widehat{y}`$, $`\widehat{z}`$ of the susceptibility tensor, since in this case M and B are collinear. The $`\widehat{y}`$ component of the torque operator Ty can be easily obtained for an applied magnetic field of the form

$$𝐁=B(\mathrm{cos}\theta \widehat{z}+\mathrm{sin}\theta \widehat{x}),$$
(14)
in Eq. 5, and is given by,

$$<𝐓_y>=g\mu _BB(<S_x>\mathrm{cos}\theta <S_z>\mathrm{sin}\theta )$$
(15)
where $`<S_\alpha >`$ = $`_{i=1}^N`$ $`<S_{i,\alpha }>`$ is the ground state expectation value of the component $`\alpha `$ of the total spin operator. We have computed $`<𝐓_y>`$ for Fe12 on the basis of the eigenvectors of the Hamiltonian in Eq. 5. When $`\theta `$ is $`\frac{\pi }{4}`$ then torque is maximal and $`<𝐓_y>`$ is just the difference between the $`\widehat{x}`$ and $`\widehat{z}`$ spin components. In Fig. 8 we have shown the variation of $`<𝐓_y>`$ with magnetic field. We can clearly observe the step behavior of the torque component which is a manifestation of the level crossing of singlet,triplet and quintet states.

## IV Summary and Outlook

We have implemented a general and efficient procedure that allows us to block factorize the spin Hamiltonian matrix based on its invariance under cyclic symmetry and parity operation. This method can be used in general for systems of other symmetries also. We have obtained the low-lying eigenvalue spectrum of ferric wheels up to $`\mathrm{Fe}_{12}`$ using the Cn rotational symmetry of the molecules. We have also analyzed the dispersion spectrum and structure factor. To reproduce the low temperature properties of ferric wheels we need to know the low-lying eigenvalue spectrum of these systems. We have also studied the dynamics of ferric wheel evolving a proper initial state. We observe nontrivial oscillation of magnetization in presence of an alternating magnetic field. We have also obtained the torque of the ferric wheels in the presence of a nonaxial magnetic field and find that the torque also exhibits step like behavior with field. Evidently a study including the effect of nonzero temperature on this oscillation is a challenging problem of future research.

Acknowledgement: We thank CSIR, India and DAE-BRNS, India for financial support.

Figure Captions

1. Representative $`\mathrm{M}_\mathrm{s}=0`$ state in (a) 6 spin-$`\frac{1}{2}`$ cluster (b) $`\mathrm{Fe}_{12}`$ wheel with all the sites having spin $`\mathrm{S}=\frac{5}{2}`$. Numbers in paranthesis correspond to the $`\mathrm{M}_\mathrm{s}`$ value at the site. Bit representation as well as the integer value are given just below the diagrams.

2. Plot of ground state energy (in unit of J) vs. inverse system size for ferric wheels up to Fe12.

3. Plot of lowest energy (in unit of J) in every total spin sector vs. the total spin (Stot) in case of Fe12. It clearly obeys Taft’s conjecture.

4. Plot of energy (in unit of J) vs. momentum vector $`k`$ for $`\mathrm{Fe}_6`$, $`\mathrm{Fe}_8`$, $`\mathrm{Fe}_{10}`$ and $`\mathrm{Fe}_{12}`$.

5. (a) Plot of static structure factor (Sq (in arbitary units)) for $`\mathrm{Fe}_6`$ in all the subspaces. Here ground state lies in B+ subspace. In (b), (c) and (d) we have plotted the static structure factor in A and B subspaces for $`\mathrm{Fe}_8`$, $`\mathrm{Fe}_{10}`$ and $`\mathrm{Fe}_{12}`$ respectively. In the doubly degenerate subspaces the weight of S(q) is much smaller than that in the ground state.

6. Plot of evolution of magnetization in presence of an alternating axial magnetic field of three different amplitudes (H0(in unit of J/$`\mathrm{}`$)).

7. Field amplitude (H0) dependence of frequency $`\mathrm{\Omega }`$ (defined in the text). Dotted line shows the exponential nature of $`\mathrm{\Omega }_{max}`$.

8. Plot of $`\widehat{y}`$ component of torque (in unit of $`\frac{J}{\mathrm{}rad}`$)with change in magnetic field applied at $`\theta =\pi /4`$.