text:
# Time-resolved Spectroscopy of the M15 X-ray Binary AC211/X2127+119

## 1 Introduction

Low Mass X-ray Binaries (LMXBs) are systems which harbour an accretion disc through which mass from a low mass secondary (M$`{}_{2}{}^{}1`$M) is accreted onto a neutron star or black hole primary. It is commonly accepted that LMXBs in the galactic disc and halo are the end product of binary evolution (see the review by Verbunt & van den Heuvel 1995). From early X-ray studies it is known that the rate of occurrence of persistently bright (L$`{}_{X}{}^{}>10^{36}`$ erg s-1) LMXBs in globular clusters is $`100`$ times higher than in other regions of our Galaxy (Katz 1975; Clark 1975). Therefore a different mechanism which enhances LMXB formation must operate in globular clusters. According to binary evolutionary models, it is currently believed that the high stellar densities and enhanced rates of star interaction in globular clusters provide a favourable environment for the formation of LMXBs via tidal capture through close encounters between neutron stars and ordinary stars (Fabian, Pringle & Rees 1975). So far, 13 persistently bright LMXBs have been discovered in globular clusters. Twelve of them, with a collection of their properties, are listed in tables 1 of Grindlay (1993) and Deutsch et al. (2000). The 13th object is the recently discovered M15 X-2 (White & Angelini 2001). Most of the globular cluster LMXBs are type-I X-ray burst sources, indicating that the compact object is an accreting neutron star. The increased effort to identify their optical counterparts is a critical part in the process of obtaining the dynamical properties of these sources (see e.g. Deutsch et al. 2000). Furthermore, their study provides important information about the evolution of the host globular cluster since close binaries can have a dramatic effect on the cluster dynamics (Elson et al. 1987; Hut et al. 1992).

X2127+119 was one of the first X-ray sources identified in a globular cluster (Giacconi et al. 1974) and the first ever globular cluster LMXB optically identified (Aurière, Le Févre & Terzan 1984; Aurière et al. 1986, Charles, Jones & Naylor 1986). Located within 2 arcsec of the crowded core of M15 (NGC 7078), the optical counterpart, AC211, is a $`V15`$ magnitude variable star. This makes AC211 one of the optically brightest LMXBs. Soon after its discovery, a periodicity of 8.5 hr was observed in the optical (Ilovaisky et al. 1987, Naylor et al. 1987) and X-rays (Hertz 1987), and initially associated with the orbital period of the system. However, Ilovaisky et al. (1993) showed through further optical photometry that the true orbital period of AC211 is 17.1 hr. The best ephemeris to date is provided by Homer & Charles (1998). The optical light curve of AC211 exhibits the largest amplitude orbital modulation of any X-ray binary (1.8 mag in U), showing a symmetrical partial primary eclipse and a broad dipping structure centred at phase 0.5 (Ilovaisky et al. 1993; Homer & Charles 1998). AC211 has been also studied in the UV (see e.g. Naylor et al. 1988, 1992) and EUV (Callanan, Drake & Fruscione 1999).

The low X-ray to optical luminosity ratio ($`L{}_{X}{}^{}/L_{opt}`$$`20`$) suggests that AC211 is an Accretion Disc Corona (ADC) source (see White & Holt 1982 for an outline of this model). The eclipse and dipping components of the light curves together with the orbital modulation of the X-ray absorption column density (Hert & Grindlay 1983, Ioannou et al. 2002a) can be easily explained if the inclination of AC211 is high. Then the secondary star and the outer asymmetric rim of the disc will partially obscure the ADC X-ray source causing the eclipse and dips, respectively. However, the detection by Ginga of an X-ray burst (with radius expansion) associated with AC211 implied that the compact object was a directly visible neutron star (Dotani et al. 1990; van Paradijs et al. 1990), a result apparently inconsistent with the ADC model. Recently, White & Angelini (2001), through observations with the Chandra X-ray Observatory have resolved X2127+119 into two X-ray sources: the LMXB associated with AC211 and a new persistent LMXB, M15 X-2. Both of these are separated by 2.7 arcsec. Since M15 X-2 is 1.6 times brighter in X-rays than AC211, its X-ray flux significantly contributes to the flux of the unresolved sources. This artificially increases the X-ray luminosity of AC211, contaminates the X-ray spectrum and reduces the amplitude of the orbital modulation in the X-ray light curve obtained in previous observations. Furthermore, since M15 X-2 shows no clear modulation in its X-ray light curve, it is reasonable to assume that it is a low-inclination LMXB and we can see the compact object directly. Thus the bright X-ray bursts observed by Ginga (Dotani et al. 1990) and RXTE (Smale 2001) can be associated with a neutron star in M15 X-2. These arguments provide a natural explanation for some of the AC211’s unusual multi-wavelength properties and further strengthen the ADC scenario.

However, despite these extensive studies, many of the fundamental physical parameters of AC211 remain unknown. This is largely due to the fact that in AC211, as in most of the persistently bright LMXBs, the optical flux generated in the disc by X-ray reprocessing overwhelms the absorption line spectrum of the secondary star. This makes it difficult to measure the radial velocity curve of the secondary. The background contamination from the surrounding field stars in M15 represents an annoying extra problem. Optical spectroscopy of AC211 by Naylor et al. (1988) showed the presence of emission at Heii $`\lambda 4686`$, Balmer lines (overpowered by the background) and Hei absorption features. From an analysis of the Hei absorption features and assuming an origin in the outer accretion disc, they claimed a systemic velocity for AC211 blue-shifted by $`150\pm 10`$ km s-1 with respect to the globular cluster. Such a high $`\gamma `$ velocity suggested tidal interaction as the mechanism of formation of the binary system accompanied by subsequent cluster ejection. However, observations by Ilovaisky et al. (1989) clearly displayed variations in the velocity of the Hei feature which showed that it was not representative of the true AC211 systemic velocity. Two alternative explanations for the Hei $`\lambda 4471`$ velocity have been suggested: Fabian, Guilbert & Callanan (1987) pointed to outflowing mass above the accretion disc as cause of the Hei absorption. In favour of this hypothesis was the discovery from HST observations of P-Cygni profiles in the H$`\beta `$ and Civ $`\lambda 1550`$ emission lines (Downes, Anderson & Margon 1996; Ioannou et al. 2002b). On the other hand, Bailyn, Garcia & Grindlay (1989) argued that the blueshifted Hei can originate in an outflowing gas from the $`L_2`$ Lagrangian point owing to a common envelope in the binary system. However, recent re-analysis of previously published optical spectroscopy by van Zyl et al. (2002) shows an orbital modulation for the Hei $`\lambda 4471`$ inconsistent with this model (see also Callanan, Naylor & Charles 1990).

Furthermore, no definitive picture of the accretion disc structure for AC211 has emerged. So far, for example, the origin of the dipping behaviour observed at phase $`0.5`$ is unknown: it may be due to a bulge caused by the re-impact of a gas stream which overflows the initial impact with the accretion disc. The lack of reliable dynamical properties for this system and observational details of the disc structure motivated us to carry out phase-resolved spectroscopy of the Heii $`\lambda 4686`$ feature and the Hei absorption lines. We proceed by giving in Section 2 a detailed discussion of the observations and data reduction. In Section 3 the analysis of the spectra is presented. A discussion of our results and a possible model that could account for some of the properties of AC211 is given in Section 4. Finally, a summary of our conclusions is presented in Section 5.

## 2 Observations and Data Reduction

AC211 was observed with the 4.5-m Multiple Mirror Telescope (MMT) on Mount Hopkins in two runs during 1996 August 10-14 UT and August 31, September 1-2 UT. Despite the variable observing conditions and mediocre seeing throughout all the observing nights (with a seeing FWHM of about 2 arcsec), the nights were good enough for spectroscopic observations. For both runs, we used the Blue Channel Spectrograph equipped with the Loral 3072x1024 CCD and the 832 grooves mm-1 grating in second order. A CuSO4 filter was used to block the first order light. This setting provided a spectral coverage of 1000 Å, a dispersion scale of 0.35 Å pix-1 and a spectral resolution of 1 Å FWHM. The slit was rotated at position angles which reduced contamination from nearby cluster stars and the spatial binning of the chip was carried out as a function of the seeing conditions. The exposure times were typically 600 and 900 seconds. HeNeAr arc lamp spectra were taken regularly through the nights in order to determine the wavelength calibration. A record of the observations and instrumental setup is provided in Table 1.

The images were bias and flat-field corrected with standard iraf111iraf is distributed by the National Optical Astronomy Observatories. routines. Spectra were extracted with the iraf kpnoslit package. The extraction of AC211 spectra and four nearby field stars was a process necessarily interactive and it was repeated several times image by image until an optimum result was obtained. The apertures were edited for each star and sized to maximize the number of pixels contributing to the object. Next, the spectra profiles were traced and the spectra were extracted without background subtraction because the slit did not include any sky region. Note that we expect a negligible sky contribution since the spectral interval observed is free of strong line bands (see e.g. Massey & Foltz 2000) and the main source of background is the unresolved globular cluster light and that from AC207 (Auriére & Cordoni 1981): the AC211 spectra presented here suffer considerably from contamination by the latter.

The pixel-to-wavelength calibration was derived from cubic spline fits to about 19 arc lines, giving a root-mean square deviation of typically $`0.03`$ Å. The wavelength scales of neighbouring arc spectra were interpolated in time and the accuracy of the whole process was checked by cross-correlation of the normalized sum of two field star spectra in the range $`\lambda \lambda 49005000`$ with the individual spectra of the corresponding star. The systematic velocity errors obtained from the scatter of the measured zero velocities were less than 3 and 5 km s-1 for the first and second run, respectively. Two spectra were rejected because they gave residual velocities of 7 and 4-$`\sigma `$.

We also found that in a few cases the trace was lost when extracting the spectra. In general, this happened near the ends of the slit where the throughput of the spectrograph drops. We removed these regions in the spectra by masking them. On the night of September 1, we suffered difficulties because of condensation inside the CCD control box which introduced spurious negative features in the object and arc frame spectra. Because of that, we rejected this night from the analysis. Our final data set consists of 40 and 19 spectra of AC211 for the first and second run respectively.

For the subsequent spectra analyses we used the MOLLY package. The spectra were rebinned into a constant velocity scale and normalized by dividing them by the result of fitting a low order spline to the continuum after masking the Balmer absorption lines and the Heii $`\lambda 4686`$ emission line in the case of AC211. To fold the spectra on the orbital period, we used the linear ephemerides from Homer et al. (1998; 2002) to determine the orbital phases:

$$\varphi =0.0,JD=2447790.964(2)\pm n\times 0.7130207(5)d$$

where $`\varphi =0.0`$ corresponds to closest approach of the secondary to the observer.

## 3 Results

### 3.1 The mean spectrum

A comparison between the mean spectrum of AC211 over the two runs is shown in Fig. 1. The only feature in the spectrum with unquestionable origin from AC211 is the Heii $`\lambda 4686`$ emission line. For the first run, this line shows a full width zero intensity (FWZI) of 865 km s-1, a FWHM=516 km s-1 and a mean equivalent width (EW) of 1.0 $`\pm `$ 0.3 Å, with the error accounting (from now on) for line variability. The mean line shows a double-peaked profile with a peak-to-peak separation of $`288\pm 6`$ km s-1 and enhanced emission on the blue-shifted peak. For the second run, we measured a FWZI=800 km s-1, a FWHM=373 km s-1 and an EW=0.6$`\pm `$ 0.1 Å for a line which shows no double peak profile. Note that the EW estimates in particular are underestimated because of contamination from AC207.

The spectral range observed was chosen to cover the Hei absorption lines reported at $`\lambda \lambda 4388,4471,4713,4921`$ by Naylor et al. (1988, 1989) and Ilovaisky (1989). A meticulous search for them in the mean and individual spectra was unsuccessful. An absorption feature at $`\lambda 4469`$ (rest frame wavelength) was discerned in the AC211 spectrum, but it is also present in the field stars (see Fig. 2), indicating that the absorption in AC211 is due to contamination from the globular cluster background. We suggest that this feature is (partial or totally) due to Tiii $`\lambda 4468.5`$. Titanium lines are detected in the red giant members of M15 (see e.g. Cohen 1979; Sneden et al. 2000) and its abundance is even enhanced for Blue Horizontal Branch stars (Behr, Cohen & McCarthy 2000). We measure an equivalent width of 0.18 $`\pm `$ 0.05 Å for the AC211 spectra and 0.19 $`\pm `$ 0.03, 0.27 $`\pm `$ 0.06, 0.24 $`\pm `$ 0.03 and 0.12 $`\pm `$ 0.02 Å for the field stars (see Auriére & Cordoni 1981) AC160, AC11, AC623 and AC256 respectively.

Most of the absorption profile for the prominent H$`\beta `$ line is also due to the background contamination, but we will show below that part of this line has its origin in the binary system.

### 3.2 Phase dependent and long-term profile variations

In order to examine the variability of the line profiles with orbital phase, we folded the Heii $`\lambda 4686`$ and the H$`\beta `$ profiles into 15 orbital phase bins by using the ephemeris of Sec 2. Fig. 3 displays the results obtained for each run in the form of trailed spectrograms where two obvious features appear. In the upper-left panel, an S-wave is present in the Heii $`\lambda 4686`$ profiles. This line component reaches maximum blueshift around phase 0.75 and maximum redshift around 0.25, with intensity increasing relative to the continuum flux between orbital phases $`0.680.94`$. The modulation in intensity explains the asymmetry in intensity of the red and blue components of the double-peaked emission for the first run when averaged over the orbit. Another interesting characteristic in the behaviour of the Heii $`\lambda 4686`$ profile is the absence of clear signs of rotational disturbance during the eclipse of the accretion disc by the donor star (which ranges in orbital phases from $`0.80.15`$; Ilovaisky et al. 1993). The lack of a disturbance in the line implies that at least some of the emission is formed mainly in non-eclipsed vertical regions on or over the accretion disc.

The radial velocity curve and EWs of the Heii $`\lambda 4686`$ are plotted in Fig. 4. The radial velocities were determined by cross-correlation of each individual profile with a Gaussian template having the same width as the mean line for the corresponding run. In the case of the first run, these velocities are modulated clearly with phase, ranging from a minimum of -197 km s-1 to a maximum of $`6`$ km s-1. It is the S-wave, which moves and varies in strength across the profile, that is the cause of these variations, while the main body of the line shows no clear modulation.

The plot of the EW versus orbital phase shows that the EW is larger during the eclipse of the accretion disc by the secondary star, with higher EWs centered around phase 0.9 and lower EWs observed at phase $`0.3`$. This behaviour suggests again that a considerable fraction of the Heii emitting region is not closely confined to the orbital plane: an intrinsically variable Heii line flux combined with an eclipse by the secondary of the continuum from the accretion disc could explain the observed modulation in the EW. Nevertheless, the possibility of apparent EW variability introduced by seeing variations may also affect these measurements.

Despite the fact that a large fraction of the H$`\beta `$ line is due to the cluster background, the folded spectra in the upper-right panel of Fig. 3 show evidence for orbital motion and, therefore, for a contribution from AC211. In particular, this line develops a blue absorption wing around orbital phase 0.25. This is in contrast with the results of van Zyl et al. (2002) who found no clear dependence of the H$`\beta `$ line on the orbital cycle.

Because of poor orbital phase coverage, the observations during the second run (lower panels of Fig. 3) are of less use in trying to measure orbital or longer term variations in the line profiles. In this regard, the trailed spectra do not sample the S-wave component at Heii $`\lambda 4686`$ or the blue absorption wing at H$`\beta `$ adequately. Fortunately some information can still be extracted since the observations during both runs match at four phase intervals ($`0.00.07`$, $`0.130.2`$, $`0.40.47`$ and $`0.530.67`$). Comparison of the Heii line profiles at the same phase intervals showed no recurrence of the profile shape. This result points to changes in the gas dynamics of the binary system over a time-scale of order a few days. Since the Rossi X-ray Timing Explorer All-Sky Monitor light curve of X2127+119 shows no significant variation in the X-ray flux/hardness ratio between both epochs of observation, we think it improbable that the differences observed above between the line profiles are due to changes in the degree of X-ray irradiation.

### 3.3 Doppler tomography

In order to gain insight into the Heii $`\lambda 4686`$ emitting regions in AC211, a tomogram of this line was constructed by using the Doppler tomography technique (Marsh and Horne 1988). We applied the maximum-entropy method of building the tomograms from the profiles obtained during the first run. Note that this technique was not used with the profiles from the second run due to the poorer orbital coverage ($`<50\%`$). Also we did not attempt to combine the data (in order to improve the orbital coverage) since there is evidence of long-term variations in the line profiles (Sec 3.2).

To obtain a reliable brightness distribution, the phase-resolved spectra acquired during the eclipse of the accretion disc by the secondary star should not be included (see e.g Marsh & Horne 1990). However, because the eclipse appears to have no effect on the velocity profiles of the Heii $`\lambda 4686`$ line, we decided to build the map using the whole data set of 40 spectra. During the process of Doppler tomography, the $`\gamma `$ velocity was adjusted to the cluster mean velocity (-107.8 km s-1, Gebhardt et al. 1997). This value is in agreement with the velocities obtained by fitting a Gaussian to the average profile. We note that such a velocity for the Heii line, in agreement with the cluster velocity, argues against a high systemic velocity for AC211. Fig. 5 displays the resultant tomogram. The map shows a faint ring of emission (a signature of emission arising from a rotating accretion disc) and a region of enhanced intensity located at approximately $`V_X=30`$ km s-1, $`V_Y=160`$ km s-1. This tomogram should be considered an approximate reconstruction of the velocity dependent line brightness distribution, since it does not strictly satisfy the assumptions that underlie the Doppler tomography technique, for example, that all motion is confined to the orbital plane or that all points in the system are equally visible at all times. In this regard, the image could also be affected by the inclusion of data observed during orbital phases $`0.350.65`$, when the accretion disc may partially mask the secondary star (Ilovaisky et al. 1993).

## 4 Discussion

Our Heii $`\lambda 4686`$ EWs are lower than the 1.7 Å reported by Bailyn, Garcia & Grindlay (1989) and considerably lower than the 9.7 Å measured by Downes et al. (1996) using HST observations. Since our observations were taken with inferior seeing than those refered to above, our spectra may be affected by higher background contamination which artificially increases the continuum. We attempted to decontaminate the MMT spectra of AC211 (out of eclipse) by subtracting different amounts of local background in an attempt to match the resulting continuum to the flux distribution of a black body of Teff=15,000 K, found by Downes et al. (1996). This resulted in Heii $`\lambda 4686`$ EWs increasing only by a factor of two, still at least a factor of 5 lower than the values measured from the HST observations (note that the MMT spectrum we selected for comparison and the HST spectrum were both obtained at similar orbital phases). This result suggests that the differences observed in the measured EWs are intrinsic to AC211 and reflect long-term dramatic changes in the emitting regions. This may be related to the changes observed between the two runs in the overall averaged line profile and points to significant changes in the form of the gas flow in the system. Also, in comparison to previous observations, no Hei absorption lines are detected in the spectra (even after our attempts at decontamination), suggesting again long-term variability in the gas regions where these lines are formed. However, because of the lack of simultaneous X-ray observations, we cannot rule out the possibility that the differences in strength of the Heii $`\lambda 4686`$ line and the presence/absence of the Hei absorption lines may be caused partially or totally by changes in the flux or hardness of the X-ray flux (in contrast to the discussion in Sec 3.2).

Keeping in mind the caveats mentioned in the previous section concerning the Doppler image, we now discuss the possible origins of the Heii $`\lambda 4686`$. In contrast to the tomograms presented by van Zyl et al. (2002), it is clear that flux from the accretion disc contributes to this line. The tomogram (Fig. 5) shows a diffuse ring with low velocities ($`200`$ km s-1). This is consistent with emission from the outer parts of a high inclination accretion disc as expected from the observed eclipses in the optical light curve. A more striking question is the origin of the enhanced emission observed at $`V_X=30`$ km s-1, $`V_Y=160`$ km s-1. This bright region produces the S-wave component in the trailed spectra of the Heii line (Fig. 4), which reaches maximum velocity near phase 0.75 and minimum near phase 0.25. In this regard, our S-wave corresponds to the so-called component B observed by van Zyl et al. (2002) in their Heii $`\lambda 4686`$ trailed spectra. Since the orbital modulation of the radial velocities is that expected for emission arising from the X-ray heated face of the secondary star, we decided to overplot on the Doppler image different outlines of the Roche lobe of the secondary star in order to constrain the physical parameters of the system. Knowing the orbital period, the position and the outline of the Roche lobe were calculated firstly by varying the inclination angle of the system between 60 and 90 degrees and assuming that the primary star is a neutron star of mass $`M_1=1.4`$ M. By doing this, we were able to match (at least partially) the Roche lobe to the emission feature in the tomogram for the cases with mass ratios ranging $`0.1<q<0.4`$ and therefore for semiamplitudes of the secondary ranging between $`190<K_2<230`$ km s-1. Fig. 5 displays the predicted positions of the Roche Lobe of the secondary star and the theoretical paths of the gas stream for an averaged mass ratio of $`q=M_2/M_1=0.3`$ and a semiamplitude of the secondary of $`K_2`$=210 km s-1. Note that the masses implied for the secondary star of $`0.1<M_2<0.6`$ M are consistent with the range of masses obtained by applying the properties of a stripped giant to the secondary of AC211 (see Homer et al. 1998). Additionally, for a fixed inclination of 90 degrees we found that for $`M_13`$ M the Roche lobe did not match the emission in the tomogram, which would suggest that the primary is a neutron star. However, these conclusions are based on the assumption that the S-wave arises from the secondary star: the fact that the S-wave component can sometimes be strong during the eclipse ingress (see upper panel Fig. 3), at the time when the X-ray irradiated face of the secondary (the expected main contributor to the S-wave component flux) points away from the observer, implies that the S-wave may not always originate on the secondary star and that the dynamical constraints obtained above should be treated with caution.

In addition, another result from the Doppler map is the absence of hotspot (the stream/disc impact region) expected at $`V_X<0`$,$`V_Y>0`$ and the lack of signs of reimpact of the gas stream in the lower quadrants ($`V_Y<0`$). The latter is an unexpected result for AC211: as stated previously, one might hope that the broad dip centred at phase 0.5 in the optical light curve could be due to occultation of the ADC/accretion disc by a bulge at the rim of the accretion disc caused by the reimpact of outflowing gas (see Sec 1). This map does not appear to support this interpretation unless the hotspot and the region of reimpact are fainter or comparable in brightness to the accretion disc at Heii $`\lambda 4686`$.

Nonetheless, it is tempting to associate the unusual spectroscopic properties of AC211 with its origin via tidal capture. Immediately after formation it is likely that the secondary was misaligned with respect to the orbital angular momentum vector. If the secondary spin and the orbital angular momentum vectors are not yet aligned, the gas stream will leave the inner Lagrangian point of the donor star with an additional kinematical contribution in the direction perpendicular to the orbital plane (see e.g. Lanzafame, Belvedere & Molteni 1994) and travel a distance above the orbital plane considerably larger than the vertical height of the accretion disc. Hence, the gas stream will be visible (generating the enhanced emission in the Doppler map), travel without interaction with the outer edge of the accretion disc, and fall back to the disc, hitting its inner side. Thus the impact of the gas stream will create a bulge which may cause the dip in the light curves centred at phase 0.5 with no observable bright spot since the hot regions are not dominant at Heii $`\lambda 4686`$ and/or they are obscured by the outer disc. In this scenario, the variations observed in the intensity of the S-wave and the long-term changes in the Heii line profiles are due to long term changes in the accretion geometry attributable in turn to a precessing secondary star. Indeed, there is some observational evidence for binary stars with one component rotating asynchronously and with a misaligned spin (Glȩbocki & Stawikowski 1997).

Although this is an interesting model, the current theories of tidal evolution of close binary systems predict that the alignment between the orbital axis and the rotational axis of a secondary with a convective envelope should occur earlier than or at a similar time to the circularization of the orbit (see Hut 1981, Eggleton & Kiseleva-Eggleton 2001). In particular, the tidal circularization time-scale for AC211 is estimated to be a few hundred years (Eggleton 2002). Therefore it is highly improbable that we are observing a misaligned secondary star unless the system has been formed recently or the secondary underfills its Roche lobe (both unlikely scenarios).

## 5 Summary and conclusions

In this paper we have presented time-resolved optical spectroscopy of AC211/X2127+119. Our Doppler tomogram of the Heii $`\lambda 4686`$ line shows enhanced emission at $`V_X=30`$ km s-1, $`V_Y=160`$ km s-1 superimposed on emission arising from the accretion disc. The phasing (if not the intensity modulation) of this emission suggests an origin on the irradiated face of the secondary. This in turn implies that AC211 is a high mass ratio system ($`0.1<q<0.4`$). Furthermore, the Doppler map suggests a neutron star primary and a low mass secondary star of mass $`0.1<M_2<0.6`$ M. We observe that the Heii $`\lambda 4686`$ shows considerable variation in its line profile on a time-scale of order a few days clearly unrelated to any orbital modulation. These have also been noted by van Zyl et al. (2002). Further observations are required in order to understand these changes in the emitting gas structure. Although a misaligned secondary (as might be expected immediately after tidal capture) may provide a natural explanation of some of the phenomena observed here, it is excluded on the basis of the short circularization time-scale of the orbit.

## 6 ACKNOWLEDGEMENTS

Use of MOLLY, DOPPLER and TRAILER developed largely by T. R. Marsh is acknowledged. We are grateful to P. Garnavich and R. P. Kirshner for generously obtaining some of the data used in this paper. We thank Craig Foltz for helpful comments during the data reduction, Lee Homer for providing the AC211 linear ephemeris, A. Recio-Blanco for confirming the presence of Tiii $`\lambda 4468.5`$ in Horizontal-Branch stars in M15 and P. Eggleton for his remarks about the possible tidal evolution of AC211. We also thank the anonymous referee for useful comments.


format: text

text:
# Exact time correlation functions for 𝑁 classical Heisenberg spins in the ‘squashed’ equivalent neighbor model

## I I. Introduction

Recently there has been a growing interest in the study of the properties of magnetic molecules. Mn12 ; Cr3tetra:1 ; Cr3tetra:2 ; Fe3tetra:1 ; Fe3tetra:2  The defining characteristic of these substances is the presence of a small cluster of magnetic ions located at the center of each molecule and surrounded by a complicated structure of non-magnetic chemical ligand groups. In general, the strength of the magnetic interaction between ions located in different molecules is negligible in comparison to the strength of their intramolecular interactions. Therefore, measurements of the magnetic properties of macroscopic samples reflect the underlying magnetic interactions within a single molecule.

The list of synthesized magnetic molecules has been constantly growing, even though most of the experimental activity has been focused on the determination of the magnetic properties of a molecule containing twelve manganese ions at its core, often referred to as Mn12. The theoretical tools currently used to describe the behavior of this relatively complicated structure are still rudimentary, and are based on a single-spin phenomenological Hamiltonian. Mn12  It is important to notice that a number of molecular structures containing smaller numbers of magnetic ions have already been synthesized. For some of these structures it is possible to perform a more detailed theoretical analysis of their magnetic behavior starting from a many-spin Hamiltonian.

Among the smaller clusters are a regular tetrahedron of $`\mathrm{Cr}^{3+}`$ ions ($`S=\frac{3}{2}`$) Cr3tetra:1 ; Cr3tetra:2 , Cr4, and a squashed tetrahedron of $`\mathrm{Fe}^{3+}`$ ions ($`S=\frac{5}{2}`$) Fe3tetra:1 ; Fe3tetra:2 , Fe4. For increasing values of the spin a description in terms of classical spins is expected to capture many of the features of the system qm:to:class:1 ; qm:to:class:2 . In the present paper we provide exact expressions for the time-dependent spin-spin correlation functions for the classical Heisenberg $`N`$-spin squashed equivalent neighbor model, which is the $`N`$-spin generalization of four classical Heisenberg spins on the corners of a squashed tetrahedron. Specific numerical results for the squashed tetrahedron case, $`N=4`$, are provided. Analogous studies have recently appeared for three spins on an isosceles triangle and on a chain isosceles , for four spins on a square ring square , and for the equivalent neighbor model of $`N`$ classical spins eq\_neighbor . Quantum time-dependent correlation functions have been computed for a dimer qm:to:class:1  and for three spins on an equilateral triangle qm:to:class:2 , and for a dimer of classical and quantum spins in a constant magnetic field dimer:field . The availability of time-dependent correlation functions is necessary to analyze neutron scattering experiments.

In Section II we define the Hamiltonian system to be studied and write the corresponding partition function. In Section III we present the constraints upon the various correlation functions. In Section IV we present our analytic results for arbitrary $`N`$. We evaluate the long-time behavior of the correlation functions at infinite temperature $`T`$, and provide analytic formulae for the low-$`T`$ modes for arbitrary $`N`$. In Section V, we present numerical results at low $`T`$ for the squashed tetrahedron, $`N=4`$. Section VI contains our conclusions. A collection of intermediate steps useful to the calculations is compiled in the Appendix.

## II II. The model

We consider $`N`$ classical spins of unit magnitude, $`|𝑺_i|=1`$, interacting according to the Hamiltonian

$$H=\frac{J_2}{2}\underset{\genfrac{}{}{0pt}{}{i,j=1}{ij}}{\overset{M}{}}𝑺_i𝑺_jJ_1𝑺_N\underset{i=1}{\overset{M}{}}𝑺_i,$$
(1)
where

$$MN12.$$
(2)
Introducing the total spin $`𝑺=_{i=1}^N𝑺_i`$ and the auxiliary variable $`𝑺_{1M}=_{i=1}^M𝑺_i=𝑺𝑺_N`$, the Hamiltonian (1) can be written as

$$H=\frac{J_1}{2}𝑺^2\frac{J_2J_1}{2}𝑺_{1M}^2,$$
(3)
where we have dropped the constant energy $`(J_1+MJ_2)/2`$. The partition function can then be calculated following the technique described in square ; ciftja:squashed . Letting $`s=|𝑺|`$ and $`x=|𝑺_{1M}|`$, one obtains

$`Z`$ $`=`$ $`{\displaystyle _0^M}𝑑x𝒟_M(x){\displaystyle _{|x1|}^{x+1}}s𝑑s\mathrm{exp}(\beta H)`$ (4)
$`=`$ $`{\displaystyle \frac{e^\alpha }{\alpha }}{\displaystyle _0^M}𝑑x𝒟_M(x)\mathrm{exp}(\alpha \gamma x^2)\mathrm{sinh}(2\alpha x),`$ (5)
where $`\beta =(k_\mathrm{B}T)^1`$, $`\alpha =\beta J_1/2`$, $`\gamma =J_2/J_1`$, and $`𝒟_M(x)`$ is the classical $`M`$-spin density of states,eq\_neighbor  which we redisplayed in Eq. (Appendix) in the Appendix.

In order to compute the time-dependent correlation functions, we first solve the classical equations of motion appropriate for the Hamiltonian, Eq. (3),

$$\dot{𝑺}_{N,1M}=J_1𝑺_{N,1M}\times 𝑺$$
(6)
and $`\dot{𝑺}=0`$, so that $`𝑺`$ is a constant of the motion. Following the technique illustrated in isosceles ; square ; eq\_neighbor , we obtain

$`𝑺_{N,1M}(t)`$ $`=`$ $`C_{N,1M}\widehat{𝒔}+A_{N,1M}`$ (7)
$`\times \left[\mathrm{cos}(st^{})\widehat{𝒙}\mathrm{sin}(st^{})\widehat{𝒚}\right],`$
where $`t^{}=J_1t`$, $`\widehat{𝒔}=𝑺/s=\widehat{𝒙}\times \widehat{𝒚}`$, $`C_N=(s^2x^2+1)/(2s)`$, $`C_{1M}=(s^2+x^21)/(2s)`$, $`A_N^2=1C_N^2`$, and $`A_N=A_{1M}`$.

We must also consider the equations of motion for the $`𝑺_i(t)`$, $`i=1,2,\mathrm{},M`$. In order to calculate the time correlation functions, symmetry allows us to choose just one of them, $`i=1`$. We then write $`𝑺_{1M}=𝑺_1+𝑺_{2M}`$, and solve

$$\dot{𝑺}_{1,2M}=J_2𝑺_{1,2M}\times 𝑺+(J_1J_2)𝑺_{1,2M}\times 𝑺_N.$$
(8)
After defining $`S_{1\pm }=S_{1x}\pm iS_{1y}`$, we obtain,

$`S_{1\pm }(t)`$ $`=`$ $`{\displaystyle \frac{A_NS_{1z0}}{C_{1M}}}\mathrm{exp}(ist^{}){\displaystyle \frac{A_N\mathrm{\Delta }S_{1z0}}{2(C_{1M}x)}}`$ (9)
$`\times \mathrm{exp}\left\{i\left[s+(1\gamma )x\right]t^{}+i\varphi _0\right\}`$
$`{\displaystyle \frac{A_N\mathrm{\Delta }S_{1z0}}{2(C_{1M}\pm x)}}`$
$`\times \mathrm{exp}\left\{i\left[s(1\gamma )x\right]t^{}i\varphi _0\right\},`$
$`S_{1z}(t)`$ $`=`$ $`S_{1z0}+\mathrm{\Delta }S_{1z0}\mathrm{cos}\left[(1\gamma )xt^{}+\varphi _0\right],`$ (10)
where $`\varphi _0`$ is an arbitrary phase, and similar equations for the components of $`𝑺_{2M}`$. After combining these equations with analogous ones for the components of $`𝑺_{1M}`$, the constants appearing in Eqs. (9) and (10) must satisfy

$`S_{1z0}`$ $`=`$ $`{\displaystyle \frac{C_{1M}}{2}}\left(1+{\displaystyle \frac{1y^2}{x^2}}\right),`$ (11)
$`(\mathrm{\Delta }S_{1z0})^2`$ $`=`$ $`{\displaystyle \frac{A_{1M}^2}{x^2}}\left[1{\displaystyle \frac{(x^2y^2+1)^2}{4x^2}}\right],`$ (12)
where $`y=|𝑺_{2M}|`$.

Previously, we solved these equations for the simplest case, $`M=2`$, for which $`y=1`$ is not a variable.isosceles  In that case, the correlation functions were obtained from the double integrals over $`x`$ and $`s`$, according to the weighting factors in Eq. (4). For $`M3`$, however, $`y`$ can vary over the entire range $`0yM1`$. Hence, for the explicit evaluation of the correlation functions with $`M3`$, it is useful to rewrite the expression of the partition function (4) in terms of a triple integral over $`s`$, $`x`$, and $`y`$,

$$Z=_0^{M1}𝒟_{M1}(y)𝑑y_{|y1|}^{y+1}𝑑x_{|x1|}^{x+1}s𝑑s\mathrm{exp}(\beta H).$$
(13)
## III III. Constraints

In this section we analyze the constraints upon the time-dependent spin-spin correlation functions

$$𝒞_{ij}(t)=𝑺_i(t)𝑺_j(0),$$
(14)
where the thermal average $`\mathrm{}`$ is performed by averaging over the arbitrary phase $`\varphi _0`$ and the variables $`s`$, $`x`$, and $`y`$, with respect to the canonical ensemble defined by Eq. (13). Due to the symmetry of the molecule, only four of the $`N(N+1)/2`$ correlation functions in Eq. (14) are distinct. We write these as $`𝒞_{11}(t)`$, $`𝒞_{12}(t)`$, $`𝒞_{1N}(t)`$, and $`𝒞_{NN}(t)`$. Conservation of the total spin adds a constraint,

$`s^2`$ $`=`$ $`𝒞_{NN}(t)+M𝒞_{11}(t)+2M𝒞_{1N}(t)`$ (15)
$`+M(M1)𝒞_{12}(t).`$
Finally, by writing the multispin correlation function $`𝑺_{1M}(t)𝑺_{1M}(0)`$ in two ways, we find a second constraint between two of the correlation functions,

$$sC_N=𝒞_{NN}(t)+M𝒞_{1N}(t),$$
(16)
where the constant $`C_N`$ is given just below Eq. (7). The two remaining independent correlation functions $`𝒞_{NN}(t)`$ and $`𝒞_{11}(t)`$ must then be calculated by explicitly substituting into Eq. (14) the time dependences obtained in Section II. For $`𝒞_{NN}(t)`$, this is relatively simple, as one can just use Eq. (7) for $`𝑺_N(t)`$, which is independent of $`y`$, to evaluate it. This leads to

$$𝒞_{NN}(t)=C_N^2+A_N^2\mathrm{cos}(st^{}),$$
(17)
which can be evaluated using the simplified weighting factors present in Eq. (4). From Eq. (16), this simplification also applies for $`𝒞_{1N}(t)`$. We note that Eq. (17) differs from the expression for the autocorrelation function in the $`N`$-spin classical Heisenberg equivalent neighbor model only by the $`x`$ dependence of the Hamiltonian,eq\_neighbor  which is irrelevant as $`T\mathrm{}`$.

The challenge is to calculate $`𝒞_{11}(t)`$. It is useful to separate the expression for $`𝒞_{11}(t)`$ into the four integrals $`I_i(t)`$ ($`i=0,\mathrm{},3)`$,

$$𝒞_{11}(t)=\underset{i=0}{\overset{3}{}}I_i(t).$$
(18)
The explicit triple integral representations of the $`I_i(t)`$ valid for arbitrary $`T`$ are given in the Appendix, where it is also shown how to reduce them to double integrals.

## IV IV. Analytic Results for Arbitrary $`N`$

### IV.1 A. Infinite temperature limit

Here we present our results for the correlation functions with general $`N`$ values as $`T\mathrm{}`$. As shown in the Appendix, in the limit $`T\mathrm{}`$, the triple integrals appearing in (18) can be reduced to single integrals. For $`N=4`$, the relevant density of states appearing in Eq. (13) is $`𝒟_2(x)=\frac{1}{2}\mathrm{\Theta }(x)\mathrm{\Theta }(2x)`$, so this reduction in the number of integrals is relatively simple. As $`T\mathrm{}`$, the different couplings appearing in the Hamiltonian become irrelevant for $`𝒞_{NN}(t)`$, so that it becomes equivalent to that of the $`N`$-spin equivalent-neighbor model,eq\_neighbor

$$\underset{T\mathrm{}}{lim}𝒞_{NN}(t)=1/N+M[\delta _N+f_N(t)],$$
(19)
where $`f_N(t)(t^{})^N`$ for $`t^{}1`$. Since as $`T\mathrm{}`$, $`s^2=N`$, $`x^2=M`$, and $`y^2=M1`$, from Eqs. (16) and (19), we have

$`\underset{T\mathrm{}}{lim}𝒞_{1N}(t)`$ $`=`$ $`1/N\delta _Nf_N(t).`$ (20)
For $`𝒞_{11}(t)`$ and $`𝒞_{12}(t)`$, even as $`T\mathrm{}`$, the situation is more complicated, as the results depend crucially upon the values of $`\gamma =J_2/J_1`$. As $`t\mathrm{}`$, the time-dependent trigonometric functions in $`I_1,I_2`$, and $`I_3`$ oscillate increasingly rapidly and yield vanishing contributions to $`lim_t\mathrm{}𝒞_{11}(t)`$, as stated by the Riemann-Lebesgue lemma lemma . Therefore, for arbitrary $`N`$,

$$\underset{t\mathrm{}}{lim}𝒞_{11}^{\gamma 1}(t)=I_0=S_{1z0}^2,$$
(21)
We note that $`I_0`$ depends upon $`N`$, and is a rather messy triple integral, but that as $`T\mathrm{}`$, can be evaluated exactly, as shown in the Appendix.

At infinite temperature one obtains for $`N=4`$,

$`\underset{\genfrac{}{}{0pt}{}{t\mathrm{}}{T\mathrm{}}}{lim}𝒞_{44}(t)={\displaystyle \frac{1}{4}}+3\delta _40.436345,`$ (22)
where $`\delta _4=(11/180)+(8/45)\mathrm{ln}20.062115`$, square ; eq\_neighbor  and

$`\underset{\genfrac{}{}{0pt}{}{t\mathrm{}}{T\mathrm{}}}{lim}𝒞_{11}^{\gamma 1}(t)`$ $``$ $`0.355496,`$ (23)
the exact expression for which is given in (Appendix) in the Appendix. In Table I in the Appendix, we also list the $`T\mathrm{}`$ values of $`lim_t\mathrm{}𝒞_{11}^{\gamma 1}(t)`$ for $`3N11`$, and compare them with the $`T\mathrm{}`$ values of $`lim_t\mathrm{}𝒞_{NN}(t)`$. We note that as $`T\mathrm{}`$, for each of these $`N`$ values, $`lim_t\mathrm{}𝒞_{11}^{\gamma 1}(t)<lim_t\mathrm{}𝒞_{NN}(t)`$. As $`T\mathrm{}`$, $`lim_t\mathrm{}𝒞_{NN}(t)`$ decreases monotonically with increasing $`N`$ to $`\frac{1}{3}`$ as $`N\mathrm{}`$.eq\_neighbor  Since $`lim_t\mathrm{}𝒞_{11}^{\gamma 1}(t)`$ also decreases monotonically with increasing $`N`$, and for $`8N11`$, its value is less than $`\frac{1}{3}`$, it appears that this inequality is likely to hold for all $`N`$ values.

We now turn to the long-time asymptotic behavior of $`𝒞_{11}(t)`$ at infinite $`T`$. Following the method described in isosceles , we first define $`\delta 𝒞_{ij}(t)𝒞_{ij}(t)lim_t\mathrm{}𝒞_{ij}(t)`$. For $`\gamma =0`$, the dominant behavior of $`lim_T\mathrm{}𝒞_{11}(t)`$ is given by $`I_3(t)`$, but for $`0\gamma 1`$, it is given by $`I_2(t)`$. At long times, $`\overline{t}1`$, where $`\overline{t}=(1\gamma )t^{}`$, one can evaluate the asymptotic behavior as $`T\mathrm{}`$ exactly. By integration by parts $`M`$ times, we find,

$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{\overline{t}1}}{lim}\delta 𝒞_{11}^{\gamma 0,1}(t)`$ $``$ $`{\displaystyle \underset{p=0}{\overset{E(M/2)}{}}}{\displaystyle \frac{A_{Mp}}{(\overline{t})^M}}f(M2p)`$ (24)
$`\times \mathrm{cos}[(M2p)\overline{t}+M\pi /2],`$
where

$`f(y)`$ $`=`$ $`1+y^2{\displaystyle \frac{(y^21)^2}{4y^3}}\mathrm{ln}\left({\displaystyle \frac{y+1}{y1}}\right)^2`$ (25)
and $`A_{Mp}`$ is given in the Appendix. Although the function $`f(y)`$ is non-analytic at $`y=1`$, it can be shown that its derivatives do not contribute to the long-time asymptotic behavior. In addition, for $`t^{}1`$, one can easily obtain the asymptotic expression of $`I_3(t)`$, leading to

$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{t^{}1}}{lim}\delta 𝒞_{11}^{\gamma =0}(t)`$ $``$ $`{\displaystyle \frac{\mathrm{sin}(t^{})}{4t^{}}}{\displaystyle _0^{M1}}𝑑y𝒟_{M1}(y)`$ (26)
$`\times y^3f(y).`$
In particular, for $`N=4`$, we obtain

$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{t^{}1}}{lim}\delta 𝒞_{44}(t)`$ $``$ $`{\displaystyle \frac{3}{4(t^{})^4}}\left[{\displaystyle \frac{3}{4}}\mathrm{cos}(4t^{})\right],`$ (27)
$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{\overline{t}1}}{lim}\delta 𝒞_{11}^{\gamma 0,1}(t)`$ $``$ $`{\displaystyle \frac{1}{8(\overline{t})^3}}\left[f(1)\mathrm{sin}(\overline{t})+f(3)\mathrm{sin}(3\overline{t})\right],`$
$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{t^{}1}}{lim}\delta 𝒞_{11}^{\gamma =0}(t)`$ $``$ $`\left({\displaystyle \frac{23}{30}}{\displaystyle \frac{9}{40}}\mathrm{ln}3\right){\displaystyle \frac{\mathrm{sin}(t^{})}{t^{}}},`$ (29)
where Eq. (27) was given previously.square ; eq\_neighbor

It is interesting to compare the present results to the analogous ones obtained for the isosceles triangle of spins, $`N=3`$. isosceles  For $`N4`$, the infinite-$`T`$, long-time behavior of $`𝒞_{11}(t)`$ for $`\gamma 0,1`$ is determined by the integral $`I_2(t)`$ given by Eq. (41). For $`N=3`$, an additional contribution to the infinite-$`T`$, long-time behavior of $`𝒞_{11}(t)`$ arises from $`I_3(t)`$ given by Eq. (42). isosceles  For $`N3`$, the correlation function $`𝒞_{11}(t)`$ for $`\gamma 0,1`$ decays slower than $`𝒞_{NN}(t)`$ \[denoted $`𝒞_{22}(t)`$ in Ref. isosceles  for $`N=3`$\], approaching its long-time asymptotic value at infinite temperature as $`(t^{})^M`$. In the limiting situation $`\gamma =0`$, corresponding for $`N=3`$ to the three-spin chain (or ‘two-pronged star’) and for $`N4`$ to an $`M`$-pronged star of spins equally coupled to a central one, as $`T\mathrm{}`$ and $`t^{}1`$, the correlation function is dominated by $`I_3(t)`$. In this case, $`𝒞_{11}(t)`$ approaches its asymptotic limit much more slowly, as $`(t^{})^1`$, as shown in Eq. (26).

### IV.2 B. Low-temperature correlation functions

At any finite temperature, it is not possible to reduce the time-dependent correlation functions to a single integral representation, even for $`N=4`$. Since the time-dependence of the integrand is a simple trigonometric function, it is convenient to compute the Fourier transforms of the $`\delta 𝒞_{ij}(t)`$, quantities which are anyhow of direct experimental relevance in neutron scattering experiments. In this case, it is then possible to express the Fourier transforms in terms of a single integral representation, which then allows a precise and fast numerical integration. We limit our numerical work to the case of the squashed tetrahedron, $`N=4`$ ($`M=3`$).

We define the Fourier transform as usual as

$$\delta \stackrel{~}{𝒞}_{ij}(\omega )=\frac{|J_1|}{\pi }_{\mathrm{}}^+\mathrm{}𝑑t\mathrm{exp}(i\omega t)\delta 𝒞_{ij}(t).$$
(30)
The position of the various peaks as a function of $`\gamma `$ may be obtained analytically in the $`T\mathrm{}`$ limit through an asymptotic evaluation of the integrals, or numerically by plotting the curves at large enough values of $`|\alpha |1/T`$. In the Appendix, we have sketched the derivation of the low-temperature mode frequencies for general $`N`$, for both FM and AFM cases. For ferromagnetic couplings, we then find,

$`\mathrm{\Omega }_1(\gamma )/J_1`$ $`=`$ $`\{\begin{array}{cc}M+1\hfill & \text{for }\gamma 1/M\hfill \\ 11/\gamma \hfill & \text{for }\gamma <1/M\hfill \end{array}`$ (31)
$`\mathrm{\Omega }_2(\gamma )/J_1`$ $`=`$ $`\{\begin{array}{cc}1+M\gamma \hfill & \text{for }\gamma 1/M\hfill \\ 0\hfill & \text{for }\gamma <1/M\hfill \end{array}`$ (32)
$`\mathrm{\Omega }_3(\gamma )/J_1`$ $`=`$ $`\{\begin{array}{cc}M|1\gamma |\hfill & \text{for }1\gamma 1/M\hfill \\ 11/\gamma \hfill & \text{for }\gamma <1/M\hfill \end{array}`$ (33)
$`\mathrm{\Omega }_4(\gamma )/J_1`$ $`=`$ $`\{\begin{array}{cc}|M(2\gamma )+1|\hfill & \text{for }\gamma 1/M\hfill \\ 2\left(11/\gamma \right)\hfill & \text{for }\gamma <1/M\text{,}\hfill \end{array}`$
and for antiferromagnetic couplings, we find,

$`\mathrm{\Omega }_1(\gamma )/|J_1|`$ $`=`$ $`\{\begin{array}{cc}\left|11/\gamma \right|\hfill & \text{for }\gamma 1/M\hfill \\ M1\hfill & \text{for }\gamma <1/M\hfill \end{array}`$ (35)
$`\mathrm{\Omega }_2(\gamma )/|J_1|`$ $`=`$ $`\{\begin{array}{cc}0\hfill & \text{for }\gamma 1/M\hfill \\ 1M\gamma \hfill & \text{for }\gamma <1/M\hfill \end{array}`$ (36)
$`\mathrm{\Omega }_3(\gamma )/|J_1|`$ $`=`$ $`\{\begin{array}{cc}\left|11/\gamma \right|\hfill & \text{for }\gamma 1/M\hfill \\ M(1\gamma )\hfill & \text{for }\gamma <1/M\hfill \end{array}`$ (37)
$`\mathrm{\Omega }_4(\gamma )/|J_1|`$ $`=`$ $`\{\begin{array}{cc}2\left|11/\gamma \right|\hfill & \text{for }\gamma 1/M\hfill \\ M(2\gamma )1\hfill & \text{for }\gamma <1/M\text{.}\hfill \end{array}`$
We remark that these formulae also apply for the isosceles triangle, $`M=2`$.isosceles

## V V. Low Temperature Numerical Results for $`N=4`$

In Fig. 1, we plot the mode frequencies $`\mathrm{\Omega }_i(\gamma )`$ relative to $`|J_1|`$, for the squashed tetrahedron case $`M=3`$. The upper and lower panels correspond to the FM and AFM cases, respectively. The circle in the upper panel of Fig. 1 denotes the absence of a zero-frequency peak at all temperatures for the regular tetrahedron. We have verified these mode frequencies by numerical evaluation of the explicit integral representations of $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ and $`\delta \stackrel{~}{𝒞}_{44}(\omega )`$. For example, in Fig. 2 we show the low-$`T`$ behavior of $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$, presented as $`\mathrm{log}_{10}[\delta \stackrel{~}{𝒞}_{11}(\omega )]`$ versus $`\omega /|J_1|`$. For the FM case with $`\gamma =0.3`$ at $`\alpha =50`$ pictured in the upper panel of Fig. 2, $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ exhibits very sharp peaks at the frequencies $`\mathrm{\Omega }_i`$, where $`\mathrm{\Omega }_i/J_1=4,1.9,2.1`$, and 6.1 for $`i=1,\mathrm{},4`$, respectively. $`\delta \stackrel{~}{𝒞}_{44}(\omega )`$ has a single sharp mode at the frequency $`\mathrm{\Omega }_1`$. This figure also shows that for $`\gamma =0.6`$, the FM $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ modes are also sharp at $`\alpha =50`$, appearing at 4, 2.8, 1.2, and 5.2, respectively, and at $`\gamma =0.9`$, they appear at 4, 3.7, 0.3, and 4.3, respectively. We note that the $`\mathrm{\Omega }_4`$ mode is much weaker in intensity than the other modes at this temperature.

For the AFM case, the modes tend to be much broader, as pictured for $`\gamma =0.6`$ in the bottom panel of Fig. 2. In this case, the low-$`T`$ mode frequencies satisfy $`\mathrm{\Omega }_i/|J_1|=2/3,0,2/3,4/3`$, so that $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_3`$ are degenerate. This degeneracy is evident in the shape of the combined mode, which appears to consist of two peaks with different widths, both centered at $`\omega /|J_1|=2/3`$. In addition, $`\mathrm{\Omega }_2`$ is a central peak, which grows in intensity as $`T`$ decreases.

For the special case of the three-pronged star, $`\gamma =0`$, the leading behaviors of the low-$`T`$ modes are presented in Fig. 3. For the FM star, pictured in the upper panel of Fig. 3, we have plotted $`n^2\delta \stackrel{~}{𝒞}_{ii}(\omega )`$ versus $`|\alpha |(\omega /J_1n)/n`$ for the largest amplitude modes $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, for $`i=1,4`$. Since $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ appear at $`\omega /J_1=4,1`$, respectively, and since the $`\mathrm{\Omega }_1`$ modes present in $`\delta \stackrel{~}{𝒞}_{44}`$ and $`\delta \stackrel{~}{𝒞}_{11}`$ are weaker than the $`\mathrm{\Omega }_2`$ mode, this presentation was chosen for clarity. Each of these modes was plotted at $`\alpha =5,10,`$ and 20, demonstrating the low-$`T`$ scaling that occurs. We also note that the $`\mathrm{\Omega }_2`$ mode in $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ drops discontinuously by many orders of magnitude (and to zero as $`T0`$) at $`\omega /J_1=1`$, as indicated by the $``$ sign. This behavior is very similar to that of the FM chain, except for the difference in the frequencies involved.isosceles

The AFM three-pronged star has parameters close to those present in the squashed tetrahedron Fe4.Fe3tetra:1 ; Fe3tetra:2  The strongest low-$`T`$ modes are pictured in the bottom panel of Fig. 3, in which we plotted $`\delta \stackrel{~}{𝒞}_{ii}(\omega )`$ versus $`|\alpha |(\omega /|J_1|n)/n`$ for $`i=1,4`$, $`n=1,2`$, and $`\alpha =5,10,20`$, and -40. Since these modes are sufficiently close in magnitude, the $`\delta \stackrel{~}{𝒞}_{ii}(\omega )`$ are not scaled in this figure. The additional modes at $`\omega /|J_1|=3,5`$ are very weak, and are not shown. As for the FM case, $`\mathrm{\Omega }_2`$ drops discontinuously by orders of magnitude at $`\omega /|J_1|=1`$, vanishing as $`T0`$. In addition, in both cases, the mode shapes approach uniform functions of $`|\alpha |(\omega /|J_1|n)`$ as $`T0`$. This behavior is actually simpler than that obtained for the AFM chain, isosceles  because in that case, the $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ modes present in $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ both approach the same frequency, $`\omega /|J_1|=1`$, as $`T0`$, making it difficult to separate them.

It is interesting to compare these findings with the simpler results in the case of a perfect tetrahedron (the equivalent neighbor model with $`N=4`$. eq\_neighbor  There only one low-$`T`$ mode is present, at $`\mathrm{\Omega }/J=4`$ in the ferromagnetic case, or at $`\mathrm{\Omega }=0`$ in the antiferromagnetic case. The low-$`T`$ scaling of these single modes was shown previously.eq\_neighbor  Allowing one spin to be coupled differently induces a splitting in the spectrum of low-$`T`$ magnons, a phenomenon which was already observed in the study of the isosceles triangle of spins. isosceles

In Fig. 4, we plot the full temperature dependence of the two primary modes present for the AFM case with $`\gamma =0.05`$, which is thought to be a better approximation to the parameters present in Fe4 than in the bottom panel of Fig. 3. Fe3tetra:1 ; Fe3tetra:2  In this figure, we show the results of calculations for both $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ and $`\delta \stackrel{~}{𝒞}_{44}(\omega )`$ at $`\alpha =0,1,5`$, and -20. At infinite $`T`$, $`\alpha =0`$, $`\delta \stackrel{~}{𝒞}_{44}(\omega )`$ exhibits a broad peak with a maximum at $`\omega /|J_1|1.7`$, and $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ has substantial weight at low frequencies, a well-defined peak at $`\omega /|J_1|0.8`$, and a small peak at $`\omega /|J_1|1.05`$. As $`T`$ is lowered, the peak in $`\delta \stackrel{~}{𝒞}_{44}(\omega )`$ develops into the sharp $`\mathrm{\Omega }_1`$ mode, approaching $`\omega /|J_1|=2`$ as $`T0`$. In addition, $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ develops into the two modes $`\mathrm{\Omega }_2`$ and $`\mathrm{\Omega }_1`$ at $`\omega /|J_1|=1.05`$ and 2, respectively. The minor peaks at $`\omega /|J_1|3.15`$ and 5.15 are too weak to show up on the scale used in this figure.

Finally, in Fig. 5 we show low-$`T`$ plots at $`N=4`$ of $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$ versus $`\omega |\alpha |^{1/2}/|J_1|`$ for the special points $`\gamma =\pm 1/3`$, corresponding to the onsets of the central peak of the mode $`\mathrm{\Omega }_2`$. In both cases, curves for $`|\alpha |=160,1280`$ are shown. Remarkably, the FM and AFM cases are nearly identical, when plotted in this manner. As for the similar scalings at the endpoints of the the parameter range of the central peak for the isosceles triangle, isosceles  this scaling only applies to the frequency, without a corresponding scaling of $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$, so that the overall scaling does not correspond to a scaling of the time in $`\delta 𝒞_{11}(t)`$. However, for the isosceles triangle, the FM and AFM cases appeared to be nearly similar at temperatures that differed by a factor of about 8, whereas for the squashed tetrahedron, the temperatures are essentially identical.

## VI VI. Conclusions

We have solved for the time correlation functions of the $`N`$-spin squashed equivalent neighbor model, with one spin coupled via the classical Heisenberg exchange $`J_1`$ to the $`M=N1`$ other spins, all of which are coupled to each other via a different Heisenberg exchange $`J_2`$. Our results are qualitatively similar to those of the isosceles triangle, $`N=3`$, but show that for arbitrary $`N3`$, there are only four low-temperature modes, given by Eqs. (31)-(IV.2) and (35)-(IV.2) for ferromagnetic and antiferromagnetic signs of $`J_1`$, respectively.

At infinite $`T`$, we showed explicitly that the long-time asymptotic behavior of the autocorrelation function $`𝒞_{11}^{\gamma =0}(t)`$ on a prong of an $`M`$-pronged star approaches its asymtotic limit as $`(t^{})^1`$. We also showed that for $`3N8`$, the infinite-$`T`$, long-time asymptotic limit of $`𝒞_{NN}(t)`$ is greater than that of $`𝒞_{11}^{\gamma 1}(t)`$, and speculate that this relation is likely to hold for arbitrary $`N`$. We also showed that at infinite $`T`$, $`𝒞_{11}^{0\gamma 1}(t)`$ approaches its long-time asymptotic limit as $`(\overline{t})^M`$, one power slower than does $`𝒞_{NN}(t)`$.

We showed explicitly that these mode frequencies apply for the isosceles triangle $`(N=3)`$ and for the squashed tetrahedron ($`N=4`$). isosceles  For the particular parameter values appropriate for the single molecule magnet Fe4, with four $`S=5/2`$ Fe+3 spins on the corners of a squashed tetrahedron, we expect that this classical calculation of the Fourier transform of the time correlation functions will represent a reasonably good envelope of the $`\delta `$-functions present in the quantum mechanical treatment of this model, provided that the temperatures are not too low with respect to $`|J_1|`$. Thus, we expect the qualitative features shown in Fig. 4 and the lower panel of Fig. 3 to be observable in inelastic neutron scattering studies of single crystals of Fe4.

## Appendix

The integrals appearing for $`M3`$ in Eq. (18) are

$`I_0`$ $`=`$ $`S_{1z0}^2,`$ (39)
$`I_1(t)`$ $`=`$ $`{\displaystyle \frac{A_N^2S_{1z0}^2}{C_{1M}^2}}\mathrm{cos}(st^{}),`$ (40)
$`I_2(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\mathrm{\Delta }S_{1z0})^2\mathrm{cos}[(1\gamma )xt^{}],`$ (41)
$`I_3(t)`$ $`=`$ $`{\displaystyle \frac{1}{4}}A_N^2(\mathrm{\Delta }S_{1z0})^2({\displaystyle \frac{\mathrm{cos}\left\{[s+(1\gamma )x]t^{}\right\}}{(x+C_{1M})^2}}`$ (42)
$`+{\displaystyle \frac{\mathrm{cos}\left\{[s(1\gamma )x]t^{}\right\}}{(xC_{1M})^2}}).`$
The density of states for $`N`$ spins is given by eq\_neighbor

$`𝒟_N(x)`$ $`=`$ $`\mathrm{\Theta }(x){\displaystyle \underset{p=0}{\overset{E[(N1)/2]}{}}}\mathrm{\Theta }(N2px)`$ (43)
$`\times \mathrm{\Theta }(xN+2p+2)d_{N2p}(x),`$
$`d_{N2p}(x)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{p}{}}}{\displaystyle \frac{(1)^k(N2kx)^{N2}}{2^{N1}(N2)!}}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{k}}\right),`$
where $`E(x)`$ is the largest integer in $`x`$ and $`\mathrm{\Theta }(x)`$ is the Heaviside step function. As noted in Eq. (13), the $`I_i(t)`$ contain integrations over $`𝒟_{M1}(y)`$.

Although the $`I_i(t)`$ are explicitly triple integrals over $`x`$, $`y`$, and $`s`$, the only $`y`$ dependence of the integrand appears in the expressions for $`S_{1z0}^2`$ and $`(\mathrm{\Delta }S_{1z0})^2`$, given by Eqs. (11) and (12) plus the expressions following Eq. (7). In most of these integrals, one has to evaluate

$`I_i`$ $`=`$ $`{\displaystyle _0^{M1}}𝑑y{\displaystyle _{|y1|}^{y+1}}𝑑xg_M(x,y)f(x,t)`$ (46)
$`=`$ $`{\displaystyle _0^{M2}}𝑑xf(x,t){\displaystyle _{|x1|}^{x+1}}𝑑yg_M(x,y)`$
$`+{\displaystyle _{M2}^M}𝑑xf(x,t){\displaystyle _{|x1|}^{M1}}𝑑yg_M(x,y),`$
where $`g_M(x,y)`$ has either the form $`a(x)[1(x^2+1y^2)^2/(4x^2)]𝒟_{M1}(y)`$ or the form $`a(x)(x^2+1y^2)^2𝒟_{M1}(y)`$, and $`f(x,t)`$ involves an integral over $`s`$. In most cases, the $`y`$ integrals can be performed before the $`x`$ integrals, reducing the triple integrals to double integrals, precisely as was done for the equivalent neighbor model with $`MN`$.eq\_neighbor

We now calculate the exact infinite-time, infinite-temperature limit of the correlation function $`𝒞_{11}(t)`$ for $`\gamma 1`$ from Eq. (39). We first perform the integration over $`s`$, and then invert the order of the remaining two integrations, as outlined above. For $`N=4`$, we then find,

$`\underset{\genfrac{}{}{0pt}{}{t\mathrm{}}{T\mathrm{}}}{lim}𝒞_{11}^{\gamma 1}(t)`$ $`=`$ $`{\displaystyle \frac{29}{360}}+{\displaystyle \frac{\pi ^2}{384}}+{\displaystyle \frac{83}{360}}\mathrm{ln}2+{\displaystyle \frac{3}{40}}\mathrm{ln}3`$
$`{\displaystyle \frac{1}{96}}\left[\mathrm{Li}_2\left({\displaystyle \frac{1}{2}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{1}{3}}\right)\right]`$
$`{\displaystyle \frac{1}{192}}\left[\mathrm{ln}\left({\displaystyle \frac{2}{3}}\right)\right]^2`$
$``$ $`0.355496,`$
where $`\mathrm{Li}_2(z)`$ is the standard dilogarithm function,

$$\mathrm{Li}_2(z)=_z^0\frac{\mathrm{ln}(1t)}{t}𝑑t.$$
(48)
The exact formulae become increasingly complicated with increasing $`N`$, so in Table I, we only list the numerical values of those additional ones for $`3N11`$, along with those of the infinite $`t,T`$ limits of $`𝒞_{NN}(t)`$.

Next, we sketch our procedure for obtaining $`lim_T\mathrm{}\delta 𝒞_{11}^{0\gamma 1}(t)`$ as $`\overline{t}1`$ for arbitrary $`N`$. From Eq. (41), we first perform the integration over $`s`$, giving us a function proportional to $`f(x)`$ given by Eq. (25). To avoid the singularity at $`x=1`$, we do not invert the order of the remaining two integrals, but instead integrate with respect to $`x`$ by parts twice, leading to

$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{\overline{t}1}}{lim}\delta 𝒞_{11}^{\gamma 0,1}(t)`$ $``$ $`{\displaystyle \frac{1}{8\overline{t}^2}}{\displaystyle \underset{\sigma =\pm 1}{}}{\displaystyle _0^{M1}}y𝑑y𝒟_{M1}(y)`$ (49)
$`\times f(y+\sigma )\mathrm{cos}[(y+\sigma )\overline{t}].`$
We then integrate with respect to $`y`$ a total of $`M2`$ times, noting that all terms proportional to derivatives of $`f`$ sum to zero. We finally obtain,

$`\underset{\genfrac{}{}{0pt}{}{T\mathrm{}}{\overline{t}1}}{lim}\delta 𝒞_{11}^{\gamma 0,1}(t)`$ $``$ $`{\displaystyle \frac{1}{(\overline{t})^M}}{\displaystyle \underset{p=0}{\overset{\mathrm{E}(M/2)}{}}}A_{Mp}f(M2p)`$ (50)
$`\times \mathrm{cos}[(M2p)\overline{t}+M\pi /2],`$
$`A_{Mp}`$ $`=`$ $`{\displaystyle \frac{(1)^{p+M}}{2^{M+1}}}[(1\delta _{p,M/2})`$ (51)
$`\times (M2p1)\left({\displaystyle \genfrac{}{}{0pt}{}{M1}{p}}\right)`$
$`(1\delta _{p,0})(M2p+1)`$
$`\times \left({\displaystyle \genfrac{}{}{0pt}{}{M1}{p1}}\right)],`$
where $`f(x)`$ is given by Eq. (25). We note that $`f(0)=\frac{8}{3}`$ and $`f(1)=2`$.

We now sketch our derivations of the low-temperature mode frequencies. We first note from Eqs. (17) and (40) that the Fourier transforms of $`\delta 𝒞_{NN}(t)`$ and $`I_1(t)`$ both contain $`\delta (s\stackrel{~}{\omega })`$, where $`\stackrel{~}{\omega }=\omega /|J_1|`$. From the above discussion, each of these then can be reduced to a single integral over $`x`$,

$$K_0(\stackrel{~}{\omega })=_{|\stackrel{~}{\omega }1|}^{\mathrm{min}(\mathrm{M},\stackrel{~}{\omega }+1)}𝑑xQ_N(x,\stackrel{~}{\omega })e^{[\alpha (\gamma 1)x^2+\stackrel{~}{\omega }^2]},$$
(52)
where $`Q_N(x,\stackrel{~}{\omega })`$ is different for $`\delta \stackrel{~}{𝒞}_{NN}(\omega )`$ and the $`\mathrm{\Omega }_1`$ mode contribution to $`\delta \stackrel{~}{𝒞}_{11}(\omega )`$. In both cases it is independent of $`\alpha `$ and $`T`$, and is therefore irrelevant to the determination of the mode frequency $`\mathrm{\Omega }_1`$ in the limit $`T0`$. The integration limits arise from the condition that the $`\delta `$-function is restricted by $`|x1|sx+1`$. For the FM case, $`\alpha >0`$, we first consider the case $`\gamma <0`$. As $`\alpha \mathrm{}`$, the integral is maximized by choosing $`x`$ to have its minimum value, $`x=|\stackrel{~}{\omega }1|`$. We then maximize the resulting expression for the exponent as a function of $`\stackrel{~}{\omega }`$, which occurs at $`\stackrel{~}{\omega }=\stackrel{~}{\omega }^{}=11/\gamma `$. For $`\gamma >1`$, the minimum $`x`$ value, $`|\stackrel{~}{\omega }1|`$, is limited for large $`\stackrel{~}{\omega }`$ by $`M`$, so $`\stackrel{~}{\omega }^{}=N`$. The crossover occurs when these frequencies are equal, $`N=11/\gamma `$, or $`\gamma =1/M`$. Setting $`\stackrel{~}{\omega }^{}=\mathrm{\Omega }_1/|J_1|`$, we thus recover Eq. (31). For the AFM case as $`T0`$, $`\alpha \mathrm{}`$, we want to minimize $`(\gamma 1)x^2+\stackrel{~}{\omega }^2`$ in the exponent. For $`\gamma >1`$, this occurs at $`x=|\stackrel{~}{\omega }1|`$, and for $`\gamma <1`$, it occurs at $`x=\stackrel{~}{\omega }+1`$. In both cases, optimizing the exponent leads to $`\stackrel{~}{\omega }^{}=|11/\gamma |`$. The latter case is restricted by the limitation $`\stackrel{~}{\omega }^{}=M1`$. The crossover between these two limits occurs at $`M=|11/\gamma |`$, or $`\gamma =1/M`$. Setting $`\stackrel{~}{\omega }^{}=\mathrm{\Omega }_1/|J_1|`$, we then recover Eq. (35).

We now focus on the integral $`I_2(t)`$, Eq. (41). We first perform the $`y`$ integral as sketched above. Then, the integral over $`s`$ does not contain any time dependence, and as $`T0`$, it is dominated by the factor $`\mathrm{exp}(\alpha s^2)`$. After integration by parts, we obtain the single integral over $`x`$, which has the form

$$I_2(t)_0^M𝑑xP_N(x)\mathrm{exp}[\alpha (\gamma x^2\pm 2x)]\mathrm{cos}[(1\gamma )xt^{}],$$
(53)
where $`P_N(x)`$ is independent of $`\alpha `$, as in Eq. (4). Fourier transformation then involves the $`\delta `$-function, $`\delta (\stackrel{~}{\omega }|1\gamma |x)`$, so that the position of the mode due to $`I_2`$ is found by optimizing the expression $`\mathrm{exp}\{\alpha [\gamma \stackrel{~}{\omega }^2/(1\gamma )^2\pm 2\stackrel{~}{\omega }/|1\gamma |]\}`$. For the FM case and $`\gamma <0`$, we maximize this function with the + sign, leading to $`\stackrel{~}{\omega }^{}=11/\gamma `$. For $`\gamma >0`$, the $`\delta `$-function was restricted by $`xM`$, leading to $`\stackrel{~}{\omega }^{}=M|1\gamma |`$. These values for $`\mathrm{\Omega }_3/|J_1|=\stackrel{~}{\omega }^{}`$ are equal at $`\gamma =1/M`$. Combining, we obtain the FM $`\mathrm{\Omega }_3`$ mode frequencies, Eq. (33). For the AFM case as $`\alpha \mathrm{}`$, we choose the $``$ sign in the above exponent, and minimize $`\gamma \stackrel{~}{\omega }^2/(1\gamma )^22\stackrel{~}{\omega }`$ in the exponent. For $`\gamma >0`$, this occurs at $`\stackrel{~}{\omega }^{}=|11/\gamma |`$. For $`\gamma <0`$, the overall exponent is bounded by $`\stackrel{~}{\omega }^{}/|1\gamma |M`$. Combining, we obtain the expressions for $`\mathrm{\Omega }_3/|J_1|`$ for the AFM case, Eq. (37).

We now turn our attention to $`I_3`$. In taking the Fourier transform, there are four $`\delta `$-functions, $`\delta [\stackrel{~}{\omega }s(1\gamma )x]`$, $`\delta [\stackrel{~}{\omega }+s+(1\gamma )x]`$, $`\delta [\stackrel{~}{\omega }+s(1\gamma )x]`$, and $`\delta [\stackrel{~}{\omega }s+(1\gamma )x]`$. These $`\delta `$-functions lead after the usual reductions of the $`y`$ integrals to the following integrals, respectively,

$`K_1(\omega )`$ $`=`$ $`{\displaystyle _{\mathrm{max}[(1\stackrel{~}{\omega })/\gamma ,(\stackrel{~}{\omega }1)/(2\gamma )]}^{\mathrm{min}[M,(\stackrel{~}{\omega }+1)/(2\gamma )]}}𝑑xR_{N1}(x,\stackrel{~}{\omega })`$ (54)
$`\times f_+(x,\stackrel{~}{\omega }),`$
$`K_2(\omega )`$ $`=`$ $`\mathrm{\Theta }(\gamma 1){\displaystyle _{\mathrm{max}[0,(1+\stackrel{~}{\omega })/\gamma ,(\stackrel{~}{\omega }1)/(\gamma 2)]}^{\mathrm{min}[M,\mathrm{\Theta }(\gamma 2)(\stackrel{~}{\omega }+1)/(\gamma 2)]}}𝑑x`$ (55)
$`\times R_{N2}(x,\stackrel{~}{\omega })f_{}(x,\stackrel{~}{\omega }),`$
$`K_3(\omega )`$ $`=`$ $`\mathrm{\Theta }(1\gamma ){\displaystyle _{\mathrm{max}[0,(1+\stackrel{~}{\omega })/\gamma ,(\stackrel{~}{\omega }+1)/(2\gamma )]}^{\mathrm{min}[M,(1\stackrel{~}{\omega })/\gamma ]}}𝑑x`$ (56)
$`\times R_{N3}(x,\stackrel{~}{\omega })f_+(x,\stackrel{~}{\omega }),`$
$`K_4(\omega )`$ $`=`$ $`{\displaystyle _{\mathrm{max}[0,(\stackrel{~}{\omega }1)/\gamma ]}^{\mathrm{min}[M,(\stackrel{~}{\omega }+1)/\gamma ]}}𝑑x`$ (57)
$`\times R_{N4}(x,\stackrel{~}{\omega })f_{}(x,\stackrel{~}{\omega }),`$
$`f_\pm (x,\stackrel{~}{\omega })`$ $`=`$ $`\mathrm{exp}\{\alpha [\stackrel{~}{\omega }^2+\gamma (\gamma 1)x^2\pm 2\stackrel{~}{\omega }(\gamma 1)x]\},`$
where the $`R_{Ni}(x,\stackrel{~}{\omega })`$ are independent of $`T`$.

We first consider the AFM case of $`K_1`$, $`\alpha \mathrm{}`$. For $`\gamma >1`$, all terms in the exponent are negative, so we need to minimize the function $`\stackrel{~}{\omega }^2+\gamma (\gamma 1)x^2+2x\stackrel{~}{\omega }(\gamma 1)`$. Setting $`x=(\stackrel{~}{\omega }1)/(2\gamma )`$, and optimizing this function with respect to $`\stackrel{~}{\omega }`$, we find that its minimum occurs at $`\stackrel{~}{\omega }^{}=2(11/\gamma )`$. For $`\gamma <0`$, the last term in the function to be minimized is negative, so we take $`x=(\stackrel{~}{\omega }+1)/(2\gamma )`$ from the upper integration limit. Optimizing the function, we find $`\stackrel{~}{\omega }^{}=2(1/\gamma 1)`$, so both $`\gamma `$ regions satisfy $`\stackrel{~}{\omega }^{}=2|11/\gamma |`$. However, this is subject to the constraint on the upper integration cutoff, which is $`(\stackrel{~}{\omega }^{}+1)/(2\gamma )=M`$, or $`\stackrel{~}{\omega }^{}=2M1M\gamma `$. These values are equal at $`\gamma =1/M`$. Altogether, $`\stackrel{~}{\omega }^{}=\mathrm{\Omega }_4/|J_1|`$ for the AFM case in Eq. (IV.2). For the FM case with $`\gamma <0`$, we take $`x=(\stackrel{~}{\omega }+1)/(2\gamma )`$, optimize, and again obtain $`\stackrel{~}{\omega }^{}=2(11/\gamma )`$. The cutoff occurs when the lower limit, $`x=(\stackrel{~}{\omega }1)/(2\gamma )`$, equals $`M`$, giving $`\stackrel{~}{\omega }^{}=2M+1M\gamma `$. The crossover occurs at $`\gamma =1/M`$, as given by Eq. (IV.2) for $`\mathrm{\Omega }_4/|J_1|`$.

Next, we consider the FM case of $`K_4`$. First for $`\alpha \mathrm{}`$, $`\gamma <0`$, it is easily seen that the exponent in $`f_{}(x\stackrel{~}{\omega })`$ is positive definite. Thus, we might expect the upper limit for $`x`$ to apply. But, this is either the cutoff, $`M`$, or a negative quantity, $`(\stackrel{~}{\omega }+1)/\gamma `$. Thus, the only positive limit is the lower cutoff, $`x=(\stackrel{~}{\omega }1)/\gamma `$, which can be positive for $`\stackrel{~}{\omega }<1`$, leading to a larger exponent than obtained by setting $`x=0`$. However, $`f_{}[(\stackrel{~}{\omega }1)/\gamma ,\stackrel{~}{\omega }]=\mathrm{exp}[(\alpha /\gamma )(\stackrel{~}{\omega }^2+\gamma 1)]`$, which for $`\alpha >0`$, $`\gamma <0`$ has a maximum at $`\stackrel{~}{\omega }^{}=0`$, corresponding to a central peak. This will be the mode frequency until $`(\stackrel{~}{\omega }1)/\gamma =M`$, the upper cutoff, resulting in $`\stackrel{~}{\omega }^{}=1+M\gamma `$. The crossover occurs at $`\gamma =1/M`$. Thus, this mode reduces to $`\mathrm{\Omega }_2/|J_1|`$ as given by Eq. (32). For the AFM case for $`\gamma >0`$, we set $`x=(\omega +1)/\gamma `$, and again we find $`f_{}[(\stackrel{~}{\omega }+1)/\gamma ,\stackrel{~}{\omega }]=\mathrm{exp}[(|\alpha |/\gamma )(\stackrel{~}{\omega }^2+\gamma 1)]`$, which has a maximum at $`\stackrel{~}{\omega }^{}=0`$. This form continues until $`(1\stackrel{~}{\omega })/\gamma =M`$, which occurs at $`\stackrel{~}{\omega }^{}=1M\gamma `$. The crossover occurs at $`\gamma =1/M`$. Thus, this gives rise to the mode $`\mathrm{\Omega }_2/|J_1|`$ in Eq. (36).

We now consider the $`K_2`$ integral. This makes a very small contribution, because of the the severe limitation that it vanishes unless $`\gamma >1`$. For the AFM case, the exponent is optimized at $`x=x^{}=\stackrel{~}{\omega }/\gamma `$, and then optimizing the mode frequency with respect to $`\stackrel{~}{\omega }`$, we find that $`\stackrel{~}{\omega }^{}=0`$, so that $`K_2`$ for AFM coupling contributes to $`\mathrm{\Omega }_2/|J_1|`$. For the FM case, the maximum exponent occurs at $`x=M`$, and from the $`\delta `$-function restrictions, we see that $`K_2`$ makes a contribution to the $`\mathrm{\Omega }_4/|J_1|`$ mode.

Finally, we discuss briefly the $`K_3`$ case, for which $`\gamma <1`$. Setting $`\gamma <0`$ for the FM case, the optimum situation is obtained when $`\stackrel{~}{\omega }^{}=0`$, so that it adds to the $`\mathrm{\Omega }_2/|J_1|`$ mode. For the AFM case, the optimum $`x`$ value is $`x^{}=\stackrel{~}{\omega }/\gamma `$, and this is restricted by $`xM`$. Hence, $`K_3`$ essentially makes a contribution to the $`\mathrm{\Omega }_2`$ mode for the AMF case, as well.