text:
# Kinetic energy of He atoms in liquid 4He-3He mixtures

## I Introduction

The microscopic static and dynamic properties of liquid 4He-3He mixtures are characterized by the interplay between the Fermi (3He) and Bose (4He) statistics, the interatomic interaction, and the quantum-mechanical zero-point motionkrotscheck ; Dobbs . Moreover, the Pauli exclusion principle strongly influences the stability of the mixturekrotscheck . Dilute solutions of 3He atoms in liquid 4He form a prototype quantum liquid as an example of an interacting boson-fermion mixture. Indeed the presence of 3He affects the condensate fraction $`n_0`$, the superfluid fraction $`\rho _s/\rho _4`$ of 4He, the individual momentum distributions, n(p),and the single-atom mean kinetic energy, $`T`$, of the two isotopes. In recent years, considerable efforts have been addressed to the understanding of microscopic static and dynamical properties in helium mixtures from both the experimental and the theoretical points of view Dobbs . Experimental Deep Inelastic Neutron Scattering (DINS) results have revealed significant and interesting discrepancies between theory and experiment as far as determination of the condensate fraction $`n_0`$ in the superfluid phase, mean kinetic energy $`T_3(x)`$ of the lighter isotope, and momentum distributions are concerned az3he4he ; wang . Here, $`x`$ is the concentration of 3He in the mixture. We stress that a substantial agreement between theories exists for the values of $`n_0`$ and $`T_3(x)`$ for low concentration mixtures bor97 ; mazzanti ; boninsegni ; lee ; ceperleymix . This discrepancies have to be compared with the remarkable agreement found between theory and experiments for pure high density liquid and solid 3HeDobbs ; senesi ; casulleras ; moroni  and pure fluid and solid 4Heglyde ; glydecond ; zoppi ; ceperley , respectively.

The single-atom mean kinetic energies $`T(x)`$ reflect the localization of the two isotopes in the mixture and are influenced by the mixture concentrationaz3he4he ; wang ; bor99 ; momentum . An important conclusion of DINS measurements in the concentration range $`0.0x1.0`$ az3he4he ; wang , is that $`T_3(x)`$ is essentially independent of $`x`$, indicating a local environment of the 3He atoms in the mixtures similar to that of pure liquid 3He. Ground state simulation techniques provide an insight into the local environment of 3He and 4He in liquid He mixtures and pure He liquids, allowing the evaluation of both partial radial distribution functions and single particle mean kinetic energies. In particular, simulation results for partial radial distribution functions, $`g_{\alpha ,\beta }(r)`$, in low concentration mixtures ($`x0.1`$), show two distinct features. The radial distribution function $`g_{3,4}(r)`$ and local density profile, $`\rho _{3,4}(r)`$, are very similar to $`g_{4,4}(r)`$, $`\rho _{4,4}(r)`$ and to the pure 4He radial distribution function ceperleymix ; fabrocini . The second feature is that $`g_{3,4}(r)`$ and $`\rho _{3,4}(r)`$ are markedly different from $`g_{3,3}(r)`$, $`\rho _{3,3}(r)`$ and from pure 3He radial distribution function. These findings support a picture where the 3He atoms experience a greater localization in the mixture with respect to pure 3He, while the 4He atoms show a microscopic structure similar to pure liquid 4He. The first feature accounts for the increased $`T_3`$ values with respect to pure liquid 3He ceperleymix ; bor97 ; mazzanti ; boninsegni ; bor99 . In the case of 4He, $`T_4`$ is on the contrary similar to the pure liquid value for $`x0`$, and a decrease of $`T_4`$ with increasing concentration is found, in agreement with experimental DINS results. It has to be stressed that the similarity of $`g_{3,4}(r)`$ and $`g_{4,4}(r)`$ does not necessarily imply similar values of $`T_3`$ and $`T_4`$ in the mixture. Moreover, since the atomic density $`n`$ in the mixtures is always larger than the atomic density of pure liquid 3He, DINS results show that $`T_3`$ is also independent of $`n`$; this, again, is in contrast with the widely-assessed density-dependence of mean kinetic energy of all quantum fluids and solids. As far as the experimental values of $`T_4(x)`$ in the $`0.0x0.4`$ range are concerned, these were found to be in agreement with microscopic calculations, resulting in a decrease of the kinetic energy of the 4He atoms with increasing concentration.

These findings motivated the present measurements, which were performed over a wider concentration range, i.e. $`x=0.00,0.10,0.35,0.65,0.90,1.00`$, and with an increased statistical accuracy than previous DINS experiments az3he4he ; wang . At present DINS is the only experimental technique which allows direct access to single-particle dynamical properties, such as the momentum distribution, $`n(\stackrel{}{p})`$, and mean kinetic energy $`T`$ evans . Experimentally this is achieved by exploiting the large values of wavevector and energy transfers involved in neutron scattering with epithermal neutrons momentum . The scattering process is well described within the framework of the Impulse Approximation (IA). In the IA the dynamical structure factor $`S(\stackrel{}{q},\omega )`$, which determines the scattered intensity, is given by:

$$S_{IA}(\stackrel{}{q},\omega )=n(\stackrel{}{p})\delta \left(\omega \frac{\mathrm{}q^2}{2M}\frac{\stackrel{}{q}\stackrel{}{p}}{M}\right)𝑑\stackrel{}{p},$$
(1)
where $`M`$ is the atomic mass of the struck nucleus. The scaling properties of the scattering law can be expressed in terms of a scaling function: $`J(y,\widehat{q})=\frac{\mathrm{}q}{M}`$ $`S_{IA}(\stackrel{}{q},\omega )`$, where $`y=\frac{M}{\mathrm{}q}\left(\omega \frac{\mathrm{}q^2}{2M}\right)`$ is the West scaling variable watson ; momentum . The function $`J(y,\widehat{q})`$, often referred to as the Neutron Compton Profile (NCP) or longitudinal momentum distribution watson ; momentum , represents the probability density distribution of $`y`$, the atomic momentum component along the direction of momentum transfer $`\widehat{q}`$. In the present case the dependence on the direction of the momentum transfer $`\widehat{q}`$ will be omitted given the absence of preferred orientations in the liquid samples. The values for $`T`$ are obtained by exploiting the second moment sum rule for $`J(y)`$ glyde ; watson :

$$_{\mathrm{}}^{\mathrm{}}y^2J(y)𝑑y=\sigma _y^2=\frac{2M}{3\mathrm{}^2}T$$
(2)
where $`\sigma _y`$ is the standard deviation of $`J(y)`$. DINS spectra from a liquid 3He-4He mixture will be then composed of two distinct contributions, one corresponding to the 4He NCP, $`J(y_4)`$, and, in a different region of the DINS spectra, one to the 3He NCP, $`J(y_3)`$. The $`J(y_{3,4})`$ functions can be separately analyzed and from their lineshape properties the momentum distributions and mean kinetic energies can be determined.

## II Experiment

The DINS measurements were carried out on the VESUVIO instrument, an inverse-geometry spectrometer operating at the ISIS pulsed neutron source (Chilton, Didcot-UK)vesuvio . On this instrument the NCP spectrum is reconstructed using the filter difference technique which consists of measuring the time of flight of the neutrons scattered by the sample; the final energy is selected by a resonant foil analyzer located between sample and detectors evans . For the present experiment the $`4.908`$ eV resonance of a 197Au foil filter was chosen. The scattered neutrons were detected by 32 glass scintillators (6Li-enriched fixed-angle elements) placed over an angular range $`115^o<`$ $`2\theta <144^o`$, yielding average wavevector transfers of $`q128`$Å-1 and $`q116`$Å-1 for 3He and 4He respectively. These large values of wavevector transfers ensured that deviations from the IA, generally described in terms of the Final State Effects (FSE) glyde  were negligible and did not affect significantly the recoil peak shapes rinat ; senesi . The corresponding average energy transfers accessed were $`\mathrm{}\omega 11eV`$ and $`\mathrm{}\omega 7eV`$ for 3He and 4He respectively. For highly absorbing 3He this energy yields a favourable ratio between the absorption and scattering cross section of about $`30`$ senesi .

The experiment has been performed along the $`1.96`$ K isotherm. Known amounts of gaseous 3He and 4He were first mixed in a reservoir at $`T=293K`$. Six mixtures were prepared, for different 3He concentrations, namely $`x=1.00,0.90,0.65,0.35,0.10,0.00`$. The mixtures were then allowed to condense into the sample cell to the homogeneous liquid phase at T= $`1.96`$ K. Special attention was paid to ensure that the liquid samples were in saturated vapour pressure conditions (SVP); this was achieved by measuring the vapour pressure of the samples on the top of the sample cell by a Baratron pressure transducer. The liquid samples were contained in a square flat aluminium cell (6 cm width, 6 cm height, 0.5 cm thickness) placed in a liquid helium flow cryostat; the sample temperatures were recorded by two Ge resistance thermometers located at the upper and lower ends of the sample cell, resulting in an average temperature $`T=1.96K\pm 0.01K`$ throughout the measurements. Experimental values of vapour pressure for each mixture composition $`x`$ were found to be in agreement with SVP data in the literature sydoriak . For each composition DINS spectra were recorded for runs lasting about 24 hours. The time-of-flight spectra were and normalized to the monitor counts by using standard routines available on VESUVIO evans . In Figure 1 a typical time of flight spectrum as recorded from a single detector for the $`x=0.35`$ mixture is displayed.

From this figure, one can note that the recoil peaks from the two different atomic masses occur at well separated positions in the time of flight spectrum. This is the case for the whole set of data in the angular range explored. Due to the high values of wavevector transfer accessed, the recoil peaks can be analyzed in wavevector spaces, i. e. $`y_3=\frac{M_3}{\mathrm{}q}\left(\omega \frac{\mathrm{}q^2}{2M_3}\right)`$ and $`y_4=\frac{M_4}{\mathrm{}q}\left(\omega \frac{\mathrm{}q^2}{2M_4}\right)`$ for 3He and 4He respectively. The fixed-angle experimental resolution, $`R_l(y_{3,4})`$, where $`l`$ is the $`lth`$ fixed-angle detector element, was determined for each detector through a standard VESUVIO experimental calibration using a lead sample. The $`R_l(y_{3,4})`$ as in previous measurements on 3He senesi ; jltp  and 4He elioepl , is well described by a Voigt function, whose parameters are $`\sigma (y_3)=0.847`$ $`\AA ^1`$, $`\frac{\mathrm{\Gamma }}{2}(y_3)=1.371`$ $`\AA ^1`$ and $`\sigma (y_4)=0.839`$ $`\AA ^1`$, $`\frac{\mathrm{\Gamma }}{2}(y_4)=1.740`$ $`\AA ^1`$, where $`\sigma (y)`$ is the standard deviation of the gaussian component and $`\frac{\mathrm{\Gamma }}{2}(y)`$ is the half width at half maximum of the lorentzian component. A parallel procedure has also been set-up using Monte Carlo neutron transport simulation codes for the VESUVIO spectrometer ddnim ; jltp ; multiplo  in order to simulate the complete moderator-sample-detector neutron transport, including multiple neutron scattering and energy dependent neutron absorption. This procedure provided simulated DINS measurements, a simulated resolution function $`R_l(y),`$ which agreed with both experimental data and experimentally calibrated $`R_l(y)`$ ddnim ; jltp ; multiplo . This ensured the reliability of the current calibration procedure, and also allows the observed $`J_{3,4}(y)`$ to be described by a convolution of the longitudinal momentum distribution and the instrumental resolution function ddnim .

### II.1 Absorption and multiple scattering correction

The effects of the 3He neutron absorption on the measured scattering from the liquid mixtures have been examined in detail, using both an analytical approach as well as a deterministic Monte Carlo simulation procedure. In 3He the neutron absorption cross section is energy-dependent with the typical 1/v variation, with a value of 5333 b for 25 meV neutrons mughabghab . The incident energy range covered by the measured 4He and 3He recoil peaks were 9-15 eV and 13-20 eV respectively (using Au absorption filters) . Although the absorption cross section of 3He has a relatively smooth variation in these ranges, we analyzed in detail the effects on the measured neutron Compton profiles. An approach for the analytic correction of absorption in strongly absorbing media was first proposed by Sears sears1 . The double-differential scattering cross section is calculated by evaluating the distribution of scattered neutrons from a sample of finite size, with the quantities $`S(q,\omega )`$, the scattering function, and $`\mathrm{\Sigma }(k)`$, the total cross section per unit volume for a neutron with wavevector k, occurring as parameters of the Boltzmann equation in the neutron transport theory sears2 . In the 3He case

$$\mathrm{\Sigma }(k)=\mathrm{\Sigma }_a(k)+\mathrm{\Sigma }_s(k)\mathrm{\Sigma }_a(k)\mathrm{\Sigma }_s(k)$$
(3)
where $`\mathrm{\Sigma }_a(k)`$ and $`\mathrm{\Sigma }_s(k)`$ are the absorption and scattering cross sections per unit volume; the beam attenuation due to multiple scattering is also found to be negligible, since the ratio of double to single scattering is of the order of $`\mathrm{\Sigma }_s(k)/2\mathrm{\Sigma }(k)`$ sears2  . In the case of a slab shaped sample and backscattering geometry, as shown in Figure 2, the double-differential cross section for single scattering is given by:

$$\frac{d^2\sigma }{d\mathrm{\Omega }dE_f}=[\frac{A}{\mathrm{\Sigma }(k_i)sec(\mathrm{\Phi }_i)+\mathrm{\Sigma }(k_f)sec(\mathrm{\Phi }_f)}]n_3\sigma \frac{k_f}{4\pi k_i}S(q,\omega )$$
(4)
where A is the surface area of the sample, $`n_3`$ is the 3He atomic number density and $`\sigma `$ is the atomic scattering cross section. The wavevector dependence of $`\mathrm{\Sigma }_a(k)`$ is expressed by $`\mathrm{\Sigma }_a(k)=n_3\sigma _a(k)=n_3\frac{4\pi }{k}b_c\mathrm{"}`$ where $`b_c\mathrm{"}`$ is the imaginary part of the coherent scattering lengthsears2 ; sears3 .

Experimental time of flight spectra from different fixed-angle detectors have been transformed into $`y`$ space and described in terms of the fixed-angle neutron Compton profiles $`J_l(y)`$ using standard VESUVIO procedures; details of the transformation of time of flight spectra in fixed-angle neutron Compton profiles are accurately described in refs. evans ; baciocco ; in order to correct for the incident-wavevector dependent absorption an analytical correction has been applied to the y-transformed data:

$$J(y)_{l,corr}=J(y)_l\frac{x}{x\sigma _3+(1x)\sigma _4}[\sigma _a(k_f)+\sigma _a(k_i)\frac{k_i^2}{k_f^2\mathrm{cos}\mathrm{\Phi }_f}]$$
(5)
where $`\sigma _a(k)`$ is the wavevector dependent absorption of 3He, and $`\sigma _3`$ and $`\sigma _4`$ are the scattering cross section for 3He and 4He respectively. In the case of inverse-geometry spectrometers such as VESUVIO, the final wavevector is constant, and for a fixed angle spectrum, Eq. 5 assumes the form

$$J(y)_{l,corr}=J_l(y)[A+Bk_i]$$
(6)
where $`A`$ and $`B`$ are constants depending on the scattering angle and the concentration $`x`$; finally, due the zero-order sum rule for J(y), $`J(y)_{l,corr}`$ (the suffix $`corr`$ is then omitted in the next sections) have been normalised to unity sears1 . The concentration dependent factor was first proposed by Hilton et al. hilton  and is particularly valid under the current experimental conditions of large wavevector and energy transfers, where cross correlation between 3He and 4He cross sections are negligible. As a complementary procedure the Monte Carlo neutron transport code DINSMS multiplo  has been also employed to evaluate multiple scattering contributions and to test the analytical absorption corrections in the mixtures. In particular this code accounts for the energy dependent absorption of the mixtures employing the following expression: :

$$n[x\sigma _3+(1x)\sigma _4]e^{n[x(\sigma _3+\sigma _a)+(1x)\sigma _4]t}dt$$
(7)
which represents the probability that a neutron travelling along the $`\widehat{t}`$ direction within the slab-shaped sample will scatter between $`t`$ and $`t+dt`$. As expected, the ratio of double to single scattering intensities varied between $`1.2\%`$ for the $`x=0.1`$ mixture and $`0.1\%`$ for the $`x=0.9`$ mixture, in agreement with theoretical predictions sears1 ; sears2 . Therefore multiple scattering corrections were neglected. Moreover, by comparing the ratio of simulated data in momentum space with and without absorption contributions included, the analytical correction factor, (Eq. 5) was recovered. An example of a correction factor, as derived from simulations for a $`x=0.35`$ mixture for a scattering angle of $`135^o`$ is presented in Figure 3.

## III Data analysis and results

The experimental spectra for each mixture composition were converted to 3He and 4He momentum space. Following this procedure a single function, $`J(y)`$, averaged over all the 32 detectors was derived for each isotope and for each composition.

The NCP were analysed by simultaneously fitting the two recoil peaks appearing in $`J(y_3)`$ ($`J(y_4)`$). The component centered at $`y=0`$ in $`J(y)`$ was fitted by a model function $`M(y)`$ broadened by the instrumental resolution. The other component, centered at negative(4He) or positive (3He) $`y`$ values were fitted using a Voigt function. The model function, $`M(y)`$, used to describe the longitudinal momentum distribution was of the form:

$$M(y)=\frac{e^{\frac{y^2}{2\sigma ^2}}}{\sqrt{2\pi \sigma ^2}}[\underset{n=2l}{\overset{\mathrm{}}{}}d_nH_n(\frac{y}{\sqrt{2\sigma ^2}})]$$
(8)
where $`H_n(\frac{y}{\sqrt{2\sigma ^2}})`$ is the *n-th* Hermite polynomial, and $`\sigma `$ and $`d_n`$ are fitting parameters; this functional form was applied by Sears searspr69  for the analysis of neutron scattering from pure liquid 4He, and a generalized form, including angular dependencies, is currently used for momentum distribution spectroscopy in hydrogen containing systems on the VESUVIO spectrometer reiterprl . We used $`d_0=1`$ and $`d_2=0`$, in order that $`M(y)`$ satisfies the following sum rules:

$$_{\mathrm{}}^{\mathrm{}}M(y)𝑑y=1$$
(9)
$$_{\mathrm{}}^{\mathrm{}}M(y)y^2𝑑y=\sigma ^2=\sigma _y^2$$
(10)
In the present case orders up to $`H_4(\frac{y}{\sqrt{2\sigma ^2}})`$ were employed. The inclusion of higher order polynomials did not result in significant improvements of the fits. As an example, in Figure 4 the $`J(y_3)`$ function for three different compositions is presented, together with the fitted lineshapes. The results for the determination of the mean kinetic energies for the two isotopes and the six mixture compositions are reported in Table 1, and are shown in Figure 5, in comparison with previous experiments and several theoretical predictions.

The present results extend the range of concentration x with respect to previous DINS measurements, with improved statistical accuracy. The 4He kinetic energy is found to decrease with concentration, while the 3He kinetic energy does not depend appreciably on the concentration.

## IV Discussion and conclusions

The results of the present experiment show that the kinetic energy of 4He atoms is strongly affected by the addition of 3He; this is expected from density and quantum statistics arguments. A remarkable agreement with previous measurements and recent theories is found, with slight discrepancies with respect to finite temperature method calculations boninsegni . On the other hand the kinetic energy of the 3He atoms appears unaffected by the presence of the higher density boson fluid, which seems to promote 3He delocalization. It is illustrative to compare the $`J(y_3)`$ spectra for the $`x=0.10`$ mixture and for the pure 3He fluid, as shown in Figure 6, there is no substantial broadening of $`J(y_3)`$ on going from the pure liquid to the low-concentration mixture.

We stress that in several papers bor97 ; mazzanti ; boninsegni  the discrepancies between experiments and theories on the determinations of $`T_3`$ have been attributed to high energy tails in the dynamic structure factors, resulting in exponential-like tails in the momentum distributions, which could be masked by background noise, or are not accounted for using inadequate model functions for $`J(y_3)`$ bor97 ; mazzanti ; boninsegni . We point out, however, that since the high energy tails are due to the repulsive part of the interatomic potential, they are present in the momentum distribution of the 4He component as well. This is well illustrated by Fig. 6 of ref. bor97 , where the momentum distributions of 3He and 4He in an $`x=0.066`$ mixture have practically the same high-momentum tails. The background noise is the same for the two isotopes in the same measurement, while the spectrometer resolution narrows on going from variables $`y_4`$ to $`y_3`$ in the same spectra; therefore one would expect an increased sensitivity to the detection of the high-momentum tails in $`J(y_3)`$. To test this hypothesis we simulated a measurement of $`J(y_3)`$ for an $`x=0.066`$ mixture in the experimental configuration of the present measurements. The momentum distribution published in ref. bor97 , $`n^{(3)}(k)`$ of 3He for the $`x=0.066`$ mixture at zero temperature was converted to $`J(y_3)`$ using the general result momentum :

$$J(y)=\frac{1}{2\pi }_{|y|}^{\mathrm{}}𝑑kkn(k)$$
(11)
This function shows pronounced high-momentum contributions due to the depletion of the Fermi sphere, the Fermi wavevector being $`k_f=0.347\AA ^1`$ . The kinetic energy value, obtained from the second moment of the reconstructed $`J(y_3)`$ was $`T_3=18.2K`$. The function was used as input for the DINSMS code, which was run for a scattering angle of $`135^o`$. The input $`J(y_3)`$ and the simulated experiment are reported in Figure 7.

The simulated spectrum was fitted by the convolution of a Voigt function, representing the spectrometer resolution, and two model functions : a simple gaussian and the Gauss-Hermite expansion introduced above (Section III) with orders up to $`H_8(\frac{y}{\sqrt{2\sigma ^2}})`$. The use of a single gaussian resulted in $`T_3=16.7\pm 0.4K`$, while the use of the Gauss-Hermite expansion resulted in $`T_3=17.0\pm 0.4K`$. This indicates that the experimentally determined kinetic energies of 3He are altered by less than 2.5 K, ruling out a strong effect of the high-momentum tails on the determination of kinetic energies.

A second test was performed to compare self-consistently the datasets available from the present measurements. The pure 3He liquid data ($`x=1.00`$) were used as the calibration measurement, and $`J(y_3)`$ was modelled to have the functional form reported in ref. azjltp , corresponding to a Fermi-like momentum distribution with a discontinuity at the Fermi momentum and high-momentum exponential tails; the resulting kinetic energy was $`T_3=11.7K`$, a value close to the theoretical predictions casulleras . This function was convoluted with a Voigt function, an ”effective” resolution to be determined from the fit to the pure liquid data. The resulting Voigt function had the following parameters:$`\sigma (y_3)=1.262`$ $`\AA ^1`$, $`\frac{\mathrm{\Gamma }}{2}(y_3)=0.9173`$ $`\AA ^1`$. This ”effective” resolution was employed to fit the $`x=0.1`$ data, and the resulting kinetic energy was $`T_3=10.47\pm 0.6K`$, confirming that no increase in kinetic energy, even in low-concentration mixtures, is indicated in the DINS experiments.

This picture suggests that the local environment of the 3He atoms remains unchanged in saturated vapour pressure liquid mixtures. Within the same picture, the $`T_3(x)`$ behaviour results also in a $`T_3(n)`$ behaviour, where $`n`$ is the total atomic density, which differs radically from the widely assessed atomic-density dependence of kinetic energy of all quantum fluids and solids. Figure 8 shows the $`T_3(n)`$ behaviour of the mixtures, as compared to $`T_3(n)`$ of the pure liquid senesi ; dimeo . The atomic density in the mixture is higher than in pure liquid 3He in equilibrium, and increases up to $`n=21.87nm^3`$ for $`x=0`$, upon adding 4He kerstead . A statistically significant departure from the density-dependence observed in the pure fluid appears for $`n`$ 19 nm-3. However we stress that in this case the simultaneous changes of density and concentration prevent a thorough picture of the density and concentration dependence of $`T_3`$. Systematic studies in several (concentration, pressure, temperature) thermodynamic states not previously investigated both experimentally and theoretically are certainly needed on the experimental and theoretical sides. For example, $`T_3(n)`$ measurements on mixtures at fixed concentration (fixed concentration and increased pressure and density) would test whether the density dependence is recovered upon approaching and crossing the liquid-solid transition. However previous and present results show unambiguously a behaviour of $`T_3`$ in the mixture which is not density dependent as for other quantum fluids and is not reproduced by any simulation studies.

The interpretation of these results is far from obvious. Two independent measurements wang ; az3he4he  and the present work, have shown unequivocally that the 3He kinetic energy in the mixtures is essentially independent on concentration and density. Further, we have shown that high energy tails cannot explain this surprising result. Given the fact that a remarkable agreement between experiment and theory has been found for the pure helium liquids, we hope that our work stimulates further theoretical and experimental effort. From the experimental point of view, we can envisage that high-resolution and high statistical quality data over an extended range of temperatures, compositions and densities, can be obtained employing novel chopper and resonance filter spectrometers at pulsed neutron sources. However, the understanding of these results remains a challenge to conventional theories of liquid isotopic helium mixtures.

###### Acknowledgements.

This research has been supported by the VESUVIO project (EU Contract Number: ERBFMGECT980142). RS kindly acknowledges J. Boronat, F. Mazzanti and A. Polls for stimulating discussions. The authors acknowledge the ISIS User Support Group for the valuable technical support during experimental measurements.

## V tables

## VI Figure captions

Figure 1. Bottom: time of flight DINS spectrum from the $`x=0.35`$ mixture for a single detector at the scattering angle $`2\theta =135^o`$. The 3He and 4He signals occur at approximately 205 $`\mu s`$ and 240 $`\mu s`$ respectively; the sample cell signal is located at about 340 $`\mu s`$. Top: wavevector and energy transfer range accessed in the bottom spectrum.

Figure 2. Schematic diagram of the scattering geometry for analytical absorption corrections. S is the slab-shaped absorbing sample. $`\mathrm{\Phi }_{i,f}`$ and $`k_{i,f}`$ are the angle of incidence (scattering) with respect to the normal to the sample, and initial (final) neutron wavevectors. Note that for the present case $`\mathrm{\Phi }_i=0`$.

Figure 3. Correction factor derived from the ratio of simulated data with and without absorption contributions.

Figure 4. $`J(y_3)`$ for the $`x=0.10,0.35,`$ and $`1.00`$ mixtures (solid circles); The peak centered at negative values of momentum corresponds to the 4He peak, centered at lower recoil energy. The solid lines are the fits to the neutron Compton profiles.

Figure 5. Single-particle kinetic energies for 3He(lower panel, solid circles) and 4He (upper panel, solid circles) respectively for the six mixtures. DINS experimental results from Azuah et al. az3he4he , right-triangles. DINS experimental results from Wang et al. wang , left-triangles. Variational calculations at zero temperature for x=0.066 from Boronat et al. bor97 , solid squares. Restricted path integral Monte Carlo calculations at $`T=2K`$ from Boninsegni et al. boninsegni , open diamonds. Variational Monte Carlo calculations at zero temperature from Lee et al. lee , solid triangles. Restricted path integral Monte Carlo calculations at $`T=2K`$ from Boninsegni et al.ceperleymix , solid diamond.

Figure 6. $`J(y_3)`$ for the $`x=0.10`$ mixture, full circles and pure 3He liquid ($`x=1.00`$), open triangles.

Figure 7. Simulated $`J(y_3)`$ for a $`x=0.066`$ mixture. Upper panel, input Neutron Compton Profile derived from ref. bor97 ; lower panel, the resulting $`J(y_3)`$ after the simulation with the DINSMS code.

Figure 8. $`T_3(n)`$ for 3He-4He liquid mixtures and pure 3He liquid; present measurements on the mixtures, solid circles; DINS experimental results from Azuah et al. az3he4he , on the mixtures, right-triangles; pure liquid 3He from Dimeo et al. dimeo ,open squares; high density pure liquid 3He from Senesi et al. senesi ,open circle.


format: text

text:
# Stable Equilibrium Based on Lévy Statistics: Stochastic Collision Models Approach

## Abstract

We investigate equilibrium properties of two very different stochastic collision models: (i) the Rayleigh particle and (ii) the driven Maxwell gas. For both models the equilibrium velocity distribution is a Lévy distribution, the Maxwell distribution being a special case. We show how these models are related to fractional kinetic equations. Our work demonstrates that a stable power-law equilibrium, which is independent of details of the underlying models, is a natural generalization of Maxwell’s velocity distribution.

There is a strong analogy between the Gaussian central limit theorem (GCLT) and the relaxation to thermal equilibrium of the Boltzmann equation (BC  and Ref. therein). However the GCLT is non-unique, which may imply that standard thermal equilibrium is non-unique. The Lévy central limit theorem (LCLT) considers the case of summation of independent identically distributed random variables with an infinite variance Feller . Hence following Montroll and Shlesinger it is natural to ask if generalized equilibrium concepts based on LCLT are meaningful Montroll . Here we start answering this question using two very different types of collision models which still reveal the same type of equilibrium. We note that Lévy statistics has many physical applications review , however its possible relation to generalized forms of equilibrium statistical mechanics is an open field of research. Recently, Bobylev and Cercignani BC , investigated a non-linear Boltzmann equation with an infinite velocity variance showing that the solution exists, and obtaining certain bounds on it. In BC  the possibility of a relation between solutions of the Boltzmann equation Ernst  and LCLT was briefly pointed out.

Our goal in this Letter is to show that a new equilibrium concept naturally emerges from old stochastic collision models. Our models demonstrate that: (i) Lévy velocity distributions serve as the natural generalization of the Maxwell velocity distribution, (ii) generalized power law equilibrium can be derived from kinetic models, there is no need to postulate a specific form of power law equilibrium, and (iii) the Lévy equilibrium obtained here possesses a certain domain of attraction, is unique, and does not depend on certain details of the underlying models.

Model 1 We consider a one dimensional tracer particle with the mass $`M`$ coupled with gas particles of mass $`m`$. The tracer particle velocity is $`V_M`$. At random times the tracer particle collides with gas particles whose velocity is denoted with $`\stackrel{~}{v}_m`$. Collisions are elastic hence from conservation of momentum and energy $`V_M^+=\xi _1V_M^{}+\xi _2\stackrel{~}{v}_m`$, where $`\xi _1=\frac{1ϵ}{1+ϵ}`$, $`\xi _2=\frac{2ϵ}{1+ϵ}`$, $`ϵm/M`$ is the mass ratio and $`V_M^+`$ $`(V_M^{})`$ is the velocity of the tracer particle after (before) a collision event. The duration of the collision events is much shorter than any other time scale in the problem. The collisions occur at a uniform rate $`R`$. The probability density function (PDF) of the gas particle velocity is $`f(\stackrel{~}{v}_m)`$. This PDF does not change during the collision process, indicating that re-collisions of the gas particles and the tracer particle are neglected.

Many works considered this type of model, imposing the condition that the gas particles are distributed according to Maxwell’s law, i.e $`f(\stackrel{~}{v}_m)`$ is Maxwellian. Since we are now investigating possible generalizations of Maxwell’s law we change this strategy and assume that $`f(\stackrel{~}{v}_m)`$ is non-Maxwellian. The goal is to see when and how the tracer particle reaches a universal equilibrium, which does not depend on the detailed shape of $`f(\stackrel{~}{v}_m)`$.

We now consider the equation of motion for the tracer particle velocity PDF $`W(V_M,t)`$ with initial conditions concentrated on $`V_M(0)`$. Standard kinetic considerations yield the linear Boltzmann equation

$$\dot{W}(V_M,t)=RW(V_M,T)+R_{\mathrm{}}^{\mathrm{}}dV_M^{}_{\mathrm{}}^{\mathrm{}}d\stackrel{~}{v}_mW(V_M^{},t)f\left(\stackrel{~}{v}_m\right)\times \delta \left(V_M\xi _1V_M^{}\xi _2\stackrel{~}{v}_m\right),$$
(1)
where the delta function gives the constrain on energy and momentum conservation in collision events. Us-usual the first (second) term in Eq. (1) describes a tracer particle leaving (entering) the velocity point $`V_M`$ at time $`t`$. Eq. (1) contains convolution integrals in velocity space hence we consider now its Fourier transform (FT). Let $`\overline{W}(k,t)`$ be the FT of the velocity PDF $`W(V_M,t)`$. Using Eq. (1), the equation of motion for $`\overline{W}(k,t)`$ is a finite difference non-local equation

$$\dot{\overline{W}}(k,t)=R\overline{W}(k,t)+R\overline{W}(k\xi _1,t)\overline{f}\left(k\xi _2\right),$$
(2)
where $`\overline{f}\left(k\right)`$ is the FT of $`f(\stackrel{~}{v}_m)`$. The solution of the equation of motion Eq. (2) is obtained by iterations

$$\overline{W}(k,t)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left(Rt\right)^n\mathrm{exp}\left(Rt\right)}{n!}e^{ikV_M(0)\xi _1^n}\mathrm{\Pi }_{i=1}^n\overline{f}\left(k\xi _1^{ni}\xi _2\right).$$
(3)
This solution has a simple interpretation. The probability that the tracer particle has collided $`n`$ times with the gas particles is given according to the Poisson law $`P_n(t)=\frac{\left(Rt\right)^n}{n!}\mathrm{exp}\left(Rt\right)`$. Let $`W_n(V_M)`$ be the PDF of the tracer particle conditioned that the particle experiences $`n`$ collision events. It can be shown that the FT of $`W_n(V_M)`$ is $`\overline{W}_n(k)=e^{ikV_M(0)\xi _1^n}\mathrm{\Pi }_{i=1}^n\overline{f}\left(k\xi _1^{ni}\xi _2\right).`$ Thus Eq. (3) is a sum over the probability of having $`n`$ collision events in time interval $`(0,t)`$ times the FT of the velocity PDF after exactly $`n`$ collision event.

In the long time limit $`\overline{W}_{eq}(k)lim_t\mathrm{}\overline{W}(k,t)`$ an equilibrium is obtained from Eq. (3). We notice that when $`Rt\mathrm{}`$, $`P_n(t)=(Rt)^n\mathrm{exp}(Rt)/n!`$ is peaked in the vicinity of $`n=Rt`$ hence it is easy to see that

$$\overline{W}_{eq}\left(k\right)=\underset{n\mathrm{}}{lim}\mathrm{\Pi }_{i=1}^n\overline{f}\left(k\xi _1^{ni}\xi _2\right).$$
(4)
In what follows we investigate properties of this equilibrium. We note that similar equilibrium can be obtained also if the collision process is not described by the Poisson law, any $`P_n(t)`$ which is peaked on $`n\mathrm{}`$ when $`t\mathrm{}`$, with (nearly) zero support for finite values of $`n`$, will exhibit this behavior.

We will consider the Rayleigh weak collision limit $`ϵ0`$. This limit is important since number of collisions needed for the tracer particle to reach an equilibrium is very large. Hence in this case we may expect the emergence of a general equilibrium concept which is not sensitive to the precise details of the velocity PDF $`f(\stackrel{~}{v}_m)`$ of the gas particles. In this limit we may also expect that in a statistical sense $`\stackrel{~}{v}_mV_M`$, hence the assumption of a uniform collision rate is reasonable in this limit.

We assume that statistical properties of gas particles velocities can be characterized with an energy scale $`T`$. Since $`T`$, $`m`$ and $`\stackrel{~}{v}_m`$ are the only variables describing the gas particle we have

$$f\left(\stackrel{~}{v}_m\right)=\frac{1}{\sqrt{T/m}}q\left(\frac{\stackrel{~}{v}_m}{\sqrt{T/m}}\right).$$
(5)
We also assume that $`f(\stackrel{~}{v}_m)`$ is an even function, as expected from symmetry. The dimensionless function $`q(x)0`$ satisfies a normalization condition $`_{\mathrm{}}^{\mathrm{}}q(x)dx=1`$, otherwise it is rather general. The scaling assumption made in Eq. (5) is very natural, since the total energy of gas particles is nearly conserved.

We first consider the case where moments of $`f\left(\stackrel{~}{v}_m\right)`$ are finite. The second moment of the gas particle velocity is $`\stackrel{~}{v}_m^2=\frac{T}{m}_{\mathrm{}}^{\mathrm{}}x^2q(x)dx`$. Without loss of generality we set $`_{\mathrm{}}^{\mathrm{}}x^2q(x)dx=1`$. The scaling behavior Eq. (5) yields $`\stackrel{~}{v}_m^{2n}=\left(\frac{T}{m}\right)^nq_{2n}`$, where the moments of $`q(x)`$ are defined according to $`q_{2n}=_{\mathrm{}}^{\mathrm{}}x^{2n}q(x)dx`$. Thus the small $`k`$ expansion of the gas particle characteristic function is

$$\overline{f}\left(k\right)=1\frac{Tk^2}{2m}+q_4\left(\frac{T}{m}\right)^2\frac{k^4}{4!}+O(k^6).$$
(6)
Inserting Eq. (6) in Eq. (4) we obtain

$$\mathrm{ln}\left[\overline{W}_{eq}\left(k\right)\right]=\frac{T}{2m}g_2\left(ϵ\right)k^2+\frac{q_43}{4!}\left(\frac{T}{m}\right)^2g_4\left(ϵ\right)k^4+O(k^6),$$
(7)
where $`g_n(ϵ)=(2ϵ)^n/[(1+ϵ)^n(1ϵ)^n]`$. The interesting thing to notice is that in the limit $`ϵ0`$, the second term on the right hand side of Eq. (7) is zero, thus $`q_4`$ is an irrelevant parameter in the problem. In similar way one can show that all terms in the expansion containing $`q_{2n}`$ with $`n>1`$ vanish in the Rayleigh limit $`ϵ0`$. Thus using Eq. (7)

$$\underset{ϵ0}{lim}\mathrm{ln}\left[\overline{W}_{eq}\left(k\right)\right]=\frac{T}{2M}k^2.$$
(8)
From Eq. (8) it is easy to see that the Maxwell velocity PDF for the tracer particle $`M`$ is obtained. Thus as expected Maxwell’s equilibrium is stable in the sense that for a large class of gas particle velocity PDFs Maxwell equilibrium is obtained.

Now we assume that $`f\left(\stackrel{~}{v}_m\right)`$ has a power law behavior, i.e., $`q(x)|x|^{(1+\alpha )}`$ when $`|x|\mathrm{}`$ and $`0<\alpha <2`$. For this case the gas particle characteristic function is

$$\overline{f}\left(k\right)=$$

$$1\frac{q_\alpha }{\mathrm{\Gamma }\left(1+\alpha \right)}\left(\frac{T}{m}\right)^{\alpha /2}|k|^\alpha +\frac{q_\beta }{\mathrm{\Gamma }\left(1+\beta \right)}\left(\frac{T}{m}\right)^{\beta /2}|k|^\beta +o(|k|^\beta )$$
(9)
where $`\alpha <\beta 2\alpha `$. $`q_\alpha `$ and $`q_\beta `$ are dimensionless numbers which depend of-course on $`q(x)`$. Without loss of generality we may set $`q_\alpha =1`$. In Eq. (9) we have used the assumption that $`f(\stackrel{~}{v}_m)`$ is even.

Inserting Eq. (9) in Eq. (4) we obtain a small $`k`$ expansion $`\mathrm{ln}\left[\overline{W}_{eq}\left(k\right)\right]`$. Taking the limit $`ϵ0`$ one can show that the terms containing $`q_\beta `$ are much smaller than the leading term $`\mathrm{ln}\left[\overline{W}_{eq}\left(k\right)\right]|k|^\alpha `$. Thus $`q_\beta `$, and in a similar way higher order coefficients, become the irrelevant parameters of the problem. Thus we find that the tracer particle equilibrium characteristic function is

$$\overline{W}_{eq}\left(k\right)\mathrm{exp}\left[\frac{2^{\alpha 1}}{\alpha \mathrm{\Gamma }\left(1+\alpha \right)}\left(\frac{T}{M}\right)^{\alpha /2}\frac{|k|^\alpha }{ϵ^{1\alpha /2}}\right],$$
(10)
thus a Lévy velocity distribution for the tracer particle is obtained. For $`\alpha 2`$ the equilibrium Eq. (10) depends on $`ϵ`$, while for the Maxwell’s case $`\alpha =2`$, the equilibrium is independent of the coupling constant $`ϵ`$. Eq. (10) implies that variance of the velocity diverges when $`\alpha <2`$. For the non-Maxwellian case the velocity distribution is characterized by the scale $`(T/M)^{\alpha /2}`$ which determines the width of the velocity distribution.

The asymptotic behavior Eq. (10) is now demonstrated using numerical examples. We consider three types of gas particle velocity PDFs, for large values of $`|v_m|\mathrm{}`$ these PDFs exhibit $`f\left(\stackrel{~}{v}_m\right)|\stackrel{~}{v}_m|^{5/2}`$, namely $`\alpha =3/2`$. Case 1 $`f(\stackrel{~}{v}_m)=N_1/(1+3^{2/3}\sqrt{m/(4T)}|\stackrel{~}{v}_m|)^{5/2}`$, Case 2, $`f(\stackrel{~}{v}_m)=N_2/[1+m\stackrel{~}{v}_m^2/(0.439T)]^{5/4}`$, where $`N_1`$ and $`N_2`$ are normalization constants. Case 3, the gas particle velocity PDF is a Lévy PDF with index $`3/2`$ whose FT is $`\overline{f}\left(k\right)=\mathrm{exp}\left[\left(\frac{T}{m}\right)^{3/4}\frac{|k|^{3/2}}{\mathrm{\Gamma }\left(5/2\right)}\right]`$.

According to our theory these power law velocity PDFs, yield a Lévy equilibrium for the tracer particle when the mass ratio becomes small, Eq. (10). In Fig. 1 we show numerically exact solution of the problem for cases (1-3). These solutions, obtained using Eq. (4) for finite values of $`ϵ`$, show a good agreement between numerical results and the asymptotic theory. The Lévy equilibrium for the tracer particle is not sensitive to precise shape of the velocity distribution of the gas particle, and hence like the Maxwell distribution is stable.

We now consider a Fokker–Planck equation which describes the evolution of the tracer particle PDF $`W(V_m,t)`$ towards the Lévy equilibrium Eq. (10). The equation is of fractional order and is obtained using a small $`ϵ`$ expansion of Eq. (1) (details to be published)

$$\frac{W(V_M,t)}{t}\frac{\overline{D}}{ϵ^{1\alpha /2}}\frac{^\alpha W(V_M,t)}{|V_M|^\alpha }+\gamma \frac{}{V_M}\left[V_MW(V_M,t)\right].$$
(11)
In Eq. (11) the Riesz fractional derivative was used review , and the dissipation term is $`\gamma =2ϵR`$. A generalized Einstein relation

$$\overline{D}=\frac{2^{\alpha 1}}{\mathrm{\Gamma }\left(1+\alpha \right)}\left(\frac{T}{M}\right)^{\alpha /2}\gamma ,$$
(12)
yields the relation between the transport coefficients $`\overline{D}`$ and $`\gamma `$. When $`\alpha =2`$ the Einstein relation is recovered. Note that Metz ; Jep ; Chechkin ; Zannette  investigated related fractional processes based on a stochastic approach (e.g. Langevin Eqs. with Lévy noise). In those investigations dissipation and fluctuations were treated as though they are independent, hence the equilibrium obtained there differs from ours.

Model 2 The question remaining is wether Lévy equilibrium a general feature, which might be obtained from other collision models. Specifically, it is interesting to see if Lévy equilibrium is compatible with a non–linear Boltzmann equation approach. Since one may suspect that the Lévy behavior obtained so far is limited to linear Boltzmann models. For this aim we investigated the one dimensional driven inelastic Maxwell model (DIMM). This model was investigated extensively in recent years in the context of inelastic gases assuming finite variance boundary conditions (see details below) Santos ; BenNaimRev . Our goal is to investigate DIMM in the quasi elastic limit showing that Lévy statistics describes the equilibrium, the Maxwell-Gauss distribution is recovered in the proper limit.

First consider the inelastic Maxwell model in the absence of external driving forces $`\dot{W}(V,t)=I(V,W)`$, where the non-linear collision integral is

$$I(V,W)W(V,t)+\frac{1}{p}_{\mathrm{}}^{\mathrm{}}W(u,t)W(\frac{vqu}{p},t)du.$$
(13)
In Eq. (13) $`p=(r+1)/2`$ and $`q=1p`$ where $`r`$ is the restitution coefficient $`0<r1`$. The kinetic scheme describes a situation where momentum is conserved during collision events, while energy is conserved only when $`r=1`$. If $`r<1`$ the steady state solution of Eq. (13) is $`W_{ss}(V)=\delta (V)`$, reflecting the loss of energy during collision events. Note that for elastic collisions $`r=1`$, any initial velocity distribution is a steady state solution. This is expected (and not informative) since two identical 1D elastic particles, exchange their velocities in collision events.

Let $`\overline{W}(k,t)`$ be the FT of $`W(V,t)`$. The boundary conditions we will consider are

$$\overline{W}(k,t)1\frac{|V|^\alpha |k|^\alpha }{\mathrm{\Gamma }\left(1+\alpha \right)},$$
(14)
for small $`k`$. Since $`W(V,t)`$ is a non-negative PDF we have $`0<\alpha 2`$. Using the Boltzmann equation (13) it is easy to show that in the elastic limit $`r=1`$, $`\frac{|V|^\alpha }{t}=0`$. For the standard case $`\alpha =2`$ considered in Santos , Eq. (14) simply reflects energy conservation, i.e. $`V^2`$ is a constant of motion. For $`\alpha <2`$, $`|V|^\alpha `$ describes the width of the probability packet, which for elastic collision is a conserved quantity.

As mentioned, when $`r<1`$ the inelastic collisions will shrink any initial probability packet to be concentrated on $`V=0`$. Similar to previous work Santos  an infinitesimal heating term is added to the equation of motion, which compensates the energy loss. We will consider the boundary conditions described in Eq. (14), while Santos  considered the Gaussian case $`\alpha =2`$. To obtain behavior compatible with Eq. (14) we consider the fractional DIMM

$$\frac{W(V,t)}{t}D_\alpha \frac{^\alpha W(V,t)}{|V|^\alpha }=I(V,W).$$
(15)
It is more convenient to consider this fractional equation in Fourier space, this yields the non-linear and non-local equation

$$\dot{\overline{W}}(k,t)+\left(1+D_\alpha |k|^\alpha \right)\overline{W}(k,t)=\overline{W}(pk,t)\overline{W}(qk,t).$$
(16)
For our aim this equation gives the definition of the fractional derivative in Eq. (15). Our aim is to investigate the steady state solution of this equation in the quasi elastic limit when $`D_\alpha 0`$ and $`r1`$. This limit is taken in such a way that the boundary condition Eq. (14) is satisfied.

Using the condition $`\frac{|V|^\alpha }{t}=0`$, and Eqs. (14, 16) we obtain

$$D_\alpha =\frac{|V|^\alpha }{\mathrm{\Gamma }\left(1+\alpha \right)}\left(1p^\alpha q^\alpha \right).$$
(17)
Without loss of generality we may set now $`|V|^\alpha =1`$. For $`\alpha =1`$ we have $`D_\alpha =0`$, while for $`0<\alpha <1`$ $`D_\alpha `$ obtains negative values. The case $`\alpha =1`$ marks the transition between a finite $`(\alpha >1)`$ and infinite $`(\alpha <1)`$ first order moment of velocity $`_{\mathrm{}}^{\mathrm{}}|V|W(V,t)dV`$. For $`\alpha <1`$ no steady state is obtained, since the dissipation due to collisions is not strong enough to compensate the heating. Our results in what follows are restricted to $`1<\alpha 2`$.

Steady state solution of Eq. (16) satisfy

$$\left(1+D_\alpha |k|^\alpha \right)\overline{W}_{ss}\left(k\right)=\overline{W}_{ss}\left(pk\right)\overline{W}_{ss}\left(qk\right).$$
(18)
An iteration method is used to obtain the solution, let $`\psi (k)\mathrm{ln}\left[\overline{W}_{ss}\left(k\right)\right]`$, and using Eq. (18) we have

$$\psi \left(k\right)=\psi \left(pk\right)+\psi \left(qk\right)\mathrm{ln}\left[1+D_\alpha |k|^\alpha \right].$$
(19)
The boundary condition Eq. (14) yields $`\psi (k)|k|^\alpha /\mathrm{\Gamma }\left(1+\alpha \right)`$. The solution of Eq. (19) is obtained using the iteration rule

$$\psi _{n+1}(k)=\psi _n\left(pk\right)+\psi _n\left(qk\right)\mathrm{ln}\left[1+D_\alpha |k|^\alpha \right],$$
(20)
where $`lim_n\mathrm{}\psi _n(k)=\psi (k)`$ and the ‘initial condition’ is $`\psi _0(k)=\mathrm{ln}\left[1+D_\alpha |k|^\alpha \right]`$. Using these rules and some algebra involving series expansions, we find

$$\psi (k)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{\left(1\right)^{n+1}|k|^{\alpha n}D_\alpha ^n}{n\left(1q^{\alpha n}p^{\alpha n}\right)},$$
(21)
where the condition $`1<\alpha 2`$ was used. Inserting Eq. (17) in Eq. (21) we obtain

$$\psi (k)=\frac{|k|^\alpha }{\mathrm{\Gamma }\left(1+\alpha \right)}+\frac{\alpha (1r)}{8}\frac{|k|^{2\alpha }}{\left[\mathrm{\Gamma }\left(1+\alpha \right)\right]^2}+O(|k|^{3\alpha })$$
(22)
The interesting thing to notice, is that the second term on the right hand side of Eq. (22) vanishes when the elastic limit $`r1`$ is considered. Inserting Eq. (17) in Eq. (21) one can show that in the elastic limit $`lim_{r1}\psi (k)=|k|^\alpha /\mathrm{\Gamma }\left(1+\alpha \right)`$. Hence the steady state characteristic function is a stretched exponential

$$\underset{r1}{lim}\overline{W}_{ss}(k)=\mathrm{exp}\left[\frac{|k|^\alpha }{\mathrm{\Gamma }\left(1+\alpha \right)}\right],$$
(23)
the inverse FT of this equation yields the symmetric stable Lévy density Feller . It is rewarding to find that Maxwell and Lévy equilibrium are obtained only in the elastic limit, thus conservation of energy in the collision events is related to the Lévy–Maxwell behavior. Far from this limit results not directly related to the Gauss–Lévy central limit theorem are obtained, Eq. (21).

To conclude, we demonstrated the relation between equilibrium properties of very different types of collision models and Lévy statistics. Thus stable behavior transcends details of individual models and hence I suspect can be found in other types of collision models. To support the idea that Lévy velocity distribution might be found in other models, we note the interesting work of Min et al Min  who used numerical simulations of a long-range interacting vortex model, and showed that distribution of velocity fields are Lévy stable.


text:
# Thermodynamics of the Farey Fraction Spin Chain

## I Introduction

Phase transitions in one-dimensional systems are unusual, essentially because, as long as the interactions are of finite range and strength, any putative ordered state at finite temperature will be disrupted by thermally induced defects, and a defect in one dimension is very effective at destroying order. Despite this, there are many examples of one-dimensional systems that do exhibit a phase transition. The Farey Fraction Spin Chain (FFSC) K-O  is one such case, which has attracted interest from both physicists and mathematicians FK ; Ka-O ; Pe . (Since this work uses some methods that may be unfamiliar to the latter, we include a paragraph at the end of this section outlining our results from a mathematical viewpoint.)

One can define the FFSC as a periodic chain of sites with two possible spin states ($`A`$ or $`B`$) at each site. This model is rigorously known to exhibit a single phase transition at temperature $`\beta _c=2`$ K-O . The phase transition itself is most unusual. The low temperature state is completely ordered K-O ; C-Kn  . In the limit of a long chain, for $`\beta >\beta _c`$, the system is either all $`A`$ or all $`B`$. Therefore the free energy $`f`$ is constant and the magnetization $`m`$ (defined via the difference in the number of spins in state $`A`$ vs. those in state $`B`$) is completely saturated over this entire temperature range. Thus, even though the system has a phase transition at finite temperature, there are no thermal effects at all in the ordered state. The same thermodynamics occurs in the Knauf spin chain (KSC) K ; C-K ; K-o ; G-K , to which the FFSC is closely related.

At temperatures above the phase transition (for $`\beta <\beta _c`$), fluctuations occur, and $`f`$ decreases with $`\beta `$. Here the system is paramagnetic, since (when the external field vanishes, see below) there is no symmetry-breaking field. Thus as the temperature increases $`m`$ jumps from its saturated value in the ordered phase to zero in the high-temperature phase C-K ; C-Kn  (see Fig. 1). (The KSC behaves similarly.)

One-dimensional models with long-range ferromagnetic interactions Aiz1 ; Aiz2  are known to exhibit a discontinuity in $`m`$ at $`\beta _c`$, but in these cases the jump in $`m`$ is less than the saturation value.

The discontinuity in $`m`$ might suggest a first-order phase transition, but in our model the behavior with temperature is different. In previous work, we proved that as a function of temperature, $`f`$ exhibits a second-order transition, and the same transition occurs in the KSC and the “Farey tree” multifractal model FK .

In beginning the research reported here, our motivation was to see whether the phase transition in the FFSC, which seems to mix first- and second-order behavior, is consistent with scaling theory. Indeed, as will be made clear, it is, in the “borderline” case when a marginal variable is present. In order to see this, we extend the definition of the FFSC to include a finite external field $`h`$. We then determine the phase diagram and free energy as a function of $`\beta `$ and $`h`$, using both rigorous and renormalization group (RG) analysis.

In the following, section II defines the model. Then, in section III we prove the existence of the free energy $`f`$ with an external field, and evaluate $`f`$ for temperatures below the phase transition. In section IV we employ renormalization group arguments to find the free energy and phase diagram for temperatures above the phase transition. Section V considers a simple model that has very similar thermodynamics but is completely solvable. Section VI summarizes our results. In the Appendix we present some arguments needed to prove the existence of $`f(\beta ,h)`$ in section III.

Since our results may be of interest to mathematicians who are unfamiliar with some of the physics employed herein, we pause to include a description of them from a more mathematical point of view. Section II defines the model and the quantities of interest. More specifically, the partition function $`Z_N`$ is a two-parameter weighted sum over the (matrices defining the) Farey fractions, and the free energy $`f`$ then follows from the limiting procedure defined in (4). The main goal of our work is to find the analytic behavior of $`f`$ as a function of the real parameters $`\beta `$, the inverse temperature (so $`\beta >0`$ is implicit), and $`h`$, the external field. Regions of parameter space for which $`f`$ is analytic are (thermodynamic) phases, and the lines of singularities that separate them are phase boundaries. In section III we prove that $`f(\beta ,h)`$ exists, and compute it exactly at low temperature (for $`\beta >\beta _c`$), which constitutes part of the ordered phase. Section IV uses renormalization group methods to determine $`f`$ at high temperatures (for $`\beta `$ near $`\beta _c`$ and $`\beta <\beta _c`$). Since this method is not rigorous, from a mathematical point of view the results should be regarded as conjectures. The main conclusions are the form of the free energy in the high-temperature phase (29, 30), the equation for the phase boundary (31, 32) and the change in magnetization $`m=f/h`$ (34) and entropy $`s=\beta ^2f/\beta `$ across the phase boundary. We also find that the ordered phase, with $`f=h`$, extends to $`\beta <\beta _c`$ when $`h`$ is sufficiently large (see Fig. 2). Section IV.3 gives predictions for the behavior of $`Z_N`$ as $`N\mathrm{}`$ near the second-order point ($`\beta =\beta _c`$ and $`h=0`$). This is related to some work in number theory, but unfortunately not yet directly. Section V examines an exactly solvable model with certain similarities to the FFSC.

## II Definition of the model

The FFSC consists of a periodic chain of $`N`$ sites with two possible spin states ($`A`$ or $`B`$) at each site. The interactions are long-range, which allows a phase transition to exist in this one-dimensional system. Let the matrices

$$M_N:=\underset{i=1}{\overset{N}{}}A^{1\sigma _i}B^{\sigma _i},\sigma _i\{0,1\},$$
(1)
where $`A:=\left(\genfrac{}{}{0pt}{}{1}{1}\genfrac{}{}{0pt}{}{0}{1}\right)`$ and $`B:=\left(\genfrac{}{}{0pt}{}{1}{0}\genfrac{}{}{0pt}{}{1}{1}\right)`$ and the dependence of $`M_N`$ on $`\{\sigma _i\}`$ has been suppressed. The energy of a particular configuration with $`N`$ spins in an external field $`h`$ is given as

$$E_N:=\mathrm{ln}(T_N)+h\left(2\underset{i=1}{\overset{N}{}}\sigma _iN\right)\mathrm{with}T_N:=\mathrm{Tr}(M_N).$$
(2)
Thus our partition function is

$$Z_N(\beta ,h)=\underset{\{\sigma _i\}}{}\mathrm{Tr}(M_N)^\beta e^{\beta h\left(2_{i=1}^N\sigma _iN\right)}.$$
(3)
This definition extends the Farey fraction spin chain model to non-vanishing external field $`h`$. Given the nature of the low-temperature $`h=0`$ system, it is natural to introduce $`h`$ in this way.

The free energy is defined as

$$f(\beta ,h):=\frac{1}{\beta }\underset{N\mathrm{}}{lim}\frac{\mathrm{ln}Z_N(\beta ,h)}{N}.$$
(4)
The existence of the free energy $`f(\beta ,h)`$ follows from simple bounds using $`f(\beta ,0)`$ (see section III below).

The definition of the FFSC is somewhat unusual. The partition function is given in terms of the energy of each possible configuration, rather than via a Hamiltonian. In fact, there is no known way to express the energy exactly in terms of the spin variables K-O . Further, numerical results indicate that when one does, the Hamiltonian has all possible even interactions (and they are all ferromagnetic), so an explicit Hamiltonian representation, even if one could find it, would be exceedingly complicated.

Note that for $`h=0`$ there are two ground states with energy $`E=\mathrm{ln}2`$. The other $`2^N2`$ states have energy $`\mathrm{ln}NENc`$, where $`c`$ is a constant. Therefore the difference between the lowest excited state energy and the ground state energy diverges as $`N\mathrm{}`$.

The phase transition in this system K-O  occurs in the following way. Divide the partition function into two terms, one due to the two ground states, and the other (call it $`Z^{}`$), due to the remaining $`2^N2`$ states. The system remains in the ground states, and $`Z^{}0`$ as $`N\mathrm{}`$, until the temperature is high enough that $`Z^{}`$ diverges with $`N`$. In section V we examine a simple model that also exhibits this feature, but is completely solvable.

Our results also apply to the KSC, which has the same thermodynamics as the FFSC model at $`h=0`$ (see FK ). An external field may be included in the KSC in exactly the same way as described above for the FFSC. The “Farey tree” model of Feigenbaum et. al. F  also has the same free energy, but it is not clear how to incorporate a field $`h`$. Our finite-size results (see section IV.3) do apply when $`h=0`$, however.

## III Free energy with an external field

In this section we show rigorously that $`f(\beta ,h)`$ exists and that

$$f(\beta ,h)=|h|,$$
(5)
for $`\beta >\beta _c`$.

For $`h>0`$ it is easy to see (from (3)) that

$$2^\beta e^{\beta hN}<Z_N(\beta ,h)<Z_N(\beta ,0)e^{\beta hN}.$$
(6)
Using the definition of the free energy then gives

$$hf(\beta ,h)f(\beta ,0)h,$$
(7)
where $`f(\beta ,h)`$ is understood to be defined via (4). Now $`f(\beta ,0)`$ is rigorously known to exist K-O . In addition, we know that $`f(\beta ,0)=0`$ for $`\beta \beta _c`$ K-O , which implies (5) for $`h>0`$ ($`h<0`$ follows similarly).

To see that $`f(\beta ,h)`$ exists for the range $`0\beta <\beta _c`$ we proceed as follows (actually, our argument applies for all $`\beta 0`$). We first show that $`\left|\frac{\mathrm{log}Z_{N+1}}{N+1}\frac{\mathrm{log}Z_N}{N}\right|0`$ as $`N\mathrm{}`$. The result then follows by use of (6). Now

$$\left|\frac{N\mathrm{log}Z_{N+1}N\mathrm{log}Z_N\mathrm{log}Z_N}{N(N+1)}\right|\left|\frac{\mathrm{log}Z_{N+1}/Z_N}{N+1}\right|+\frac{1}{N+1}\left|\frac{\mathrm{log}Z_N}{N}\right|,$$

and we see by (6) and the existence of $`f(\beta ,0)`$ that the second term $`\frac{1}{N+1}\left|\frac{\mathrm{log}Z_N}{N}\right|\frac{K}{N+1}`$ for some finite constant $`K`$. In the appendix we show that $`2^\beta e^{\beta |h|}\frac{Z_{N+1}}{Z_N}2e^{\beta |h|}`$ which completes our proof of the existence of the free energy for all $`\beta 0`$ and $`h`$.

We also know rigorously that $`f(t,0)c\frac{t}{\mathrm{ln}t}+..`$, where $`c>0`$, $`t=\frac{\beta _c}{\beta }1`$, for $`t>0`$ (see Fig. 1). It follows that $`f(t,h)`$ must have at least one singularity between the regions with low and high temperatures, i.e. a phase transition from the ordered to the high-temperature phase.

Since we can not calculate $`f(\beta ,h)`$ exactly for $`\beta <\beta _c`$ (except for $`h=0`$ and $`\beta \beta _c`$ ), we use another method, in the next section, to examine the thermodynamics.

## IV  Renormalization group analysis

### IV.1 Mean field theory

In mean field theory one assumes that there is an expansion of the free energy of the form

$$f_{MF}=a+btM^2+uM^4ghM+\mathrm{},$$
(8)
where $`M`$ is the magnetization and the “constants” $`a`$, $`b`$, $`u`$ and $`g`$ are weakly dependent on the reduced temperature $`t`$ (defined at the end of section III) and external field $`h`$. Note that $`u>0`$ is required for stability, and $`b>0,g0`$ in the high-temperature phase. (The possibility that $`g=0`$ is ruled out below.)

Minimizing (8) with respect to $`M`$, one obtains the free energy and magnetization in mean field approximation. Explicitly

1. for $`t>0`$ and $`h0`$ the magnetization

$$M_0\frac{1}{6}\left[\frac{u}{gh}+\left(\frac{2bt}{3gh}\right)^3\right]^{\frac{1}{3}}$$

(note the limiting cases $`M_00`$ for $`h=0`$ and $`M_01/6(gh/u)^{1/3}`$ for $`t=0`$)
2. for $`t<0`$ and $`h0`$, but $`h`$ sufficiently small, the magnetization

$$M_0\left(\frac{b|t|}{2u}\right)^{\frac{1}{2}}+\frac{gh}{4b|t|}$$

(however when $`\left(\frac{gh}{2u}\right)^2+4\left(\frac{bt}{6u}\right)^3>0`$, $`M_0`$ is given by the $`t>0`$ formula). We include this second case only for completeness. Since our system is completely saturated at low temperatures this result is not employed in our analysis.

In the following we use the first result in an RG analysis.

### IV.2 Renormalization group analysis

We assume two relevant fields ($`t`$ and $`h`$) and one marginal field ($`u`$). These assumptions are reasonable, since our model has an Ising-like ordered state, the interactions are (apparently) all ferromagnetic, and there is a logarithmic term in the free energy.

The infinitesimal renormalization group transformation for the singular part of the free energy is

$$f_s(t,h,u)=e^d\mathrm{}f_s(t(\mathrm{}),h(\mathrm{}),u(\mathrm{})).$$
(9)
Because of the marginal field $`u`$, the analysis is somewhat more complicated than otherwise. We follow the treatment of Cardy (Ca , see also Wegner We ). The RG equations take the form

$`du/d\mathrm{}`$ $`=`$ $`xu^2+\mathrm{}`$ (10)
$`dt/dl`$ $`=`$ $`y_ttz_tut+\mathrm{}`$ (11)
$`dh/dl`$ $`=`$ $`y_hhz_huh+\mathrm{},`$ (12)
where we keep only the most important terms. The omitted terms are either higher order or go to zero more rapidly with $`\mathrm{}`$ than those included. From (10) we find (note $`t=t(0)`$, $`h=h(0)`$, $`u=u(0)`$)

$$u(\mathrm{})=\frac{u(0)}{1+xu(0)\mathrm{}}.$$
(13)
Both $`t`$ and $`h`$ have the same functional form, namely

$$\mathrm{ln}(t(\mathrm{}_0)/t(0))=y_t\mathrm{}_0\frac{z_t}{x}\mathrm{ln}[1+xu(0)\mathrm{}_0]$$
(14)
and

$$\mathrm{ln}(h(\mathrm{}_0)/h(0))=y_h\mathrm{}_0\frac{z_h}{x}\mathrm{ln}[1+xu(0)\mathrm{}_0],$$
(15)
where $`\mathrm{}_0`$ is such that $`t(\mathrm{}_0)=O(1)`$ or $`h(\mathrm{}_0)=O(1)`$. From (14) we can write

$$\mathrm{}_0\frac{1}{y_t}\mathrm{ln}\frac{t_0}{t}+\frac{z_t}{xy_t}\mathrm{ln}\left[1+\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right],$$
(16)
where we assume $`t_0/t1`$. This result together with (9) gives us

$$f_s(t,h,u)\left|\frac{t}{t_0}\right|^{\frac{d}{y_t}}\left[1+\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^{\frac{z_td}{y_tx}}f_s(t(\mathrm{}_0),h(\mathrm{}_0),u(\mathrm{}_0)).$$
(17)
Since the free energy on the rhs is evaluated at $`\mathrm{}_0`$, which is far from the critical point, it can be calculated from mean field theory. Above the critical temperature ($`t>0`$) with small external field $`h`$ ($`h(\mathrm{}_0)t(\mathrm{}_0)`$) we obtain for the free energy

$$f_s(t(\mathrm{}_0),h(\mathrm{}_0),u(\mathrm{}_0))a\frac{3(gh(\mathrm{}_0))^2}{16bt(\mathrm{}_0)}.$$
(18)
The relation between $`h(\mathrm{}_0)`$ and $`t(\mathrm{}_0)`$ follows from (14) and (15). Eliminating $`h(\mathrm{}_0)`$ allows us to rewrite (18) as

$$f_sa\left|\frac{t_0}{t}\right|^{2\frac{y_h}{y_t}}h^2\left[1+\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^{2y_h\left[\frac{z_t}{y_tx}\frac{z_h}{y_hx}\right]}\left(\frac{3g^2}{16bt(\mathrm{}_0)}\right).$$
(19)
Substituting the result into (9) gives two terms,

$$\left|\frac{t}{t_0}\right|^{\frac{d}{y_t}}\left[1+\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^{\frac{z_td}{y_tx}}a,$$
(20)
and

$$\left|\frac{t}{t_0}\right|^{\frac{d}{y_t}2\frac{y_h}{y_t}}h^2\left[1+\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^{\frac{z_td}{y_tx}+2y_h\left[\frac{z_t}{y_tx}\frac{z_h}{y_hx}\right]}\left(\frac{3g^2}{16bt(\mathrm{}_0)}\right).$$
(21)
The first term can be compared with the exact result at $`h=0`$ (see section I). It follows that

$$\frac{d}{y_t}=1=\frac{z_t}{x}.$$
(22)
The second term gives us the dependence on external field. Eliminating $`t(\mathrm{}_0)`$ instead of $`h(\mathrm{}_0)`$ we obtain

$$\frac{1}{t}\left|\frac{h}{h_0}\right|^{\frac{d}{y_h}+\frac{y_t}{y_h}}\left[1+\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^{\frac{z_hd}{y_hx}y_t\left[\frac{z_h}{y_hx}\frac{z_t}{y_tx}\right]}\left(\frac{3(gh(\mathrm{}_0))^2}{16b}\right).$$
(23)
Equating the two expressions (21) and (23) for the same term in the free energy gives us the RG eigenvalues

$$\frac{d}{y_t}=\frac{d}{y_h}=1,$$
(24)
where $`d`$ is the dimensionality of the system. This is of course one for our model, but since none of our results require setting $`d=1`$ we leave it unspecified.

Finally we can write down the singular part of the free energy for the high-temperature phase

$$f_s(t,h,u)\left|\frac{t}{t_0}\right|\left[\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^1a\frac{h^2}{t}\left[\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^{1\frac{z_h}{x}}\left(\frac{3g^2}{16b}\right).$$
(25)
Since $`f<0`$ for $`h=0`$ in this phase, (25) implies that $`a<0`$.

For the ordered phase we know rigorously that the free energy has no temperature dependence for $`h=0`$. The spins are all up or all down. When we add an external field it will break the symmetry and all the spins will be oriented in the field direction. Thus the free energy at $`\mathrm{}_0`$ is

$$f_s(t(\mathrm{}_0),h(\mathrm{}_0),u(\mathrm{}_0))=|h|(\mathrm{}_0).$$
(26)
Proceeding as in the derivation of (17) from (9) and (14) we get

$$f_s(t,h,u)\left|\frac{h}{h_0}\right|^{\frac{d}{y_h}}\left[1+\frac{x}{y_h}u\mathrm{ln}\frac{h_0}{h}\right]^{\frac{z_hd}{y_hx}}f_s(t(\mathrm{}_0),h(\mathrm{}_0),u(\mathrm{}_0)),$$
(27)
using (15), (26) and (24) then give

$$f_s(t,h,u)|h|\left[1+\frac{x}{y_h}u\mathrm{ln}\frac{h_0}{h}\right]^{\frac{z_h}{x}}.$$
(28)
Because the magnetization in the ordered state is completely saturated the logarithmic correction must vanish. Therefore $`z_h=0`$.

Thus the asymptotic form for the free energy of the high-temperature state is

$$f_s(t,h,u)\left|\frac{t}{t_0}\right|\left[\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^1a\frac{h^2}{t}\left[\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]\left(\frac{3g^2}{16b}\right).$$
(29)
We can recast this result more suggestively as

$$f_s(t,h,u)f_s(t,0,u)\frac{1}{2}h^2\chi (t,0,u),$$
(30)
where $`\chi =^2f/h^2`$ is the susceptibility. Note that $`\chi 1/f_s`$ which is consistent with scaling theory, since (using (24)), $`ft^{2\alpha }=t^{d/y_t}=t`$ while $`\chi t^\gamma =t^{(d2y_h)/y_t}=t^1`$. This relation holds regardless of whether we set the dimensionality $`d=1`$ or not. In addition, the coefficient of $`\frac{t}{\mathrm{ln}t}`$ for the free energy at $`h=0`$ and $`t0,t>0`$ is known exactly Diss ; P-S , so that the combination of constants $`\frac{y_ta}{|t_0|xu}`$ may be determined.

The phase boundary is given by the continuity of the free energy. Now we expect the ordered phase to exist for $`\beta <\beta _c`$ if $`h`$ is large enough (this is reflected in the assumption of two relevant fields-if another phase intervened there would be more). Thus one must equate the two expressions for $`f`$. One finds that the phase boundary between the ordered and high-temperature phase, close to the critical point, follows

$$|h|k\frac{t}{\mathrm{ln}t/t_0},$$
(31)
where $`k=\left\{\frac{8by_t}{3xug^2}\left[1\sqrt{1+\frac{3ag^2}{4bt_0}}\right]\right\}`$. Since $`f`$ is quadratic in $`h`$ in the high-temperature phase, there are in general two solutions with $`h>0`$. However, the one at larger $`h`$ is not physical since it gives rise to a magnetization $`m>1`$ and violates the convexity of the free energy as well, so we employ the other.

In order to find the change in magnetization across the phase boundary we use (31) with constants included

$$|h|\frac{t}{\mathrm{ln}\frac{t_0}{t}}\left\{\frac{8by_t}{3xug^2}\left[1\sqrt{1+\frac{3ag^2}{4bt_0}}\right]\right\}.$$
(32)
In arriving at (32), we (as mentioned) chose the root that makes $`m<1`$ in the high-temperature phase. Note that in the limiting case that $`\frac{3ag^2}{4bt_0}=1`$, $`m=1`$ but the two roots coincide.

Now from (29)

$$m\frac{h}{t}\left[\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]\left(\frac{3g^2}{8b}\right).$$
(33)
Eliminating the external field using (32), and since the magnetization in the ordered phase takes the values $`m\pm 1`$, we find

$$\mathrm{\Delta }m\sqrt{1+\frac{3ag^2}{4bt_0}}.$$
(34)
Note that $`t_0`$ is a constant of order one and recall that $`a<0`$, thus on the phase boundary the discontinuity in magnetization is constant (and less than one), at least close to the second-order point (we argue below that $`g=0`$ is not possible in this model). Now we can look at the change in entropy (per site) $`s=\beta ^2f/\beta `$ across the phase boundary. We get

$$\mathrm{\Delta }s2\left[\frac{x}{y_t}u\mathrm{ln}\frac{t_0}{t}\right]^1\left(\frac{a}{t_0}+\frac{4b}{3g^2}\left[1\sqrt{1+\frac{3ag^2}{4bt_0}}\right]\right).$$
(35)
These results show that the phase transition is first-order everywhere except at $`h=0`$.

In the limiting case when $`\frac{3ag^2}{4bt_0}=1`$, already mentioned, one finds that both $`\mathrm{\Delta }m=0`$ and $`\mathrm{\Delta }s=0`$. However, it is easy to see that both the susceptibility $`\chi `$ and the specific heat will have a discontinuity across the phase boundary.

Note that the magnetization change given by (34) exhibits a kind of “discontinuity of the discontinuity”, in that its limiting value as one approaches the second-order point is not the same as its value at that point. This is not the case for the entropy change, or for these quantities in the model examined in section V.

Finally, we argue that $`g=0`$ is not possible in the high-temperature phase. Since the second derivative of $`f`$ with respect to $`h`$ at $`h=0`$ is proportional to both $`g`$ and the susceptibility $`\chi `$, it suffices to demonstrate that $`\chi >0`$. It is straightforward to show that $`\chi `$ is proportional to $`\mathrm{\Sigma }_{j=1}^Ns_1s_j`$ where the spin variables $`s_i:=2\sigma _i1,s_i\{1,1\}`$ (cf. (1)), and the angular brackets denote a thermal average. Now the $`j=1`$ term in this sum is $`1`$, and due to the ferromagnetic interactions in the spin chain, the remaining terms cannot be negative. Note that this argument is not completely rigorous, since for the FFSC we only have numerical evidence that the interactions are all ferromagnetic. The KSC, on the other hand, is known to have all interactions ferromagnetic C-Kn , so that $`s_1s_j>0`$ follows from the GKS inequalities.

### IV.3 Finite-size scaling

We can use our results to make some predictions about finite-size (i.e. $`N`$ large but $`N<\mathrm{}`$) effects on the thermodynamics. We make the standard assumption that the size of our spin chain is a relevant field with eigenvalue $`1`$. Of course, since our system has long-range interactions the validity of finite-size scaling may be questioned Ca , but it is still interesting to see the results. The treatment is the same as in the case of the relevant fields $`t`$ and $`h`$. The renormalization equation for the inverse size $`I:=N^1`$ is then

$$dI/dl=Iz_IuI+\mathrm{}.$$
(36)
Thus we get

$$f_s(t,h,u,N^1)\left|\frac{N_0}{N}\right|^d\left[1+xu\mathrm{ln}\frac{N}{N_0}\right]^{\frac{z_Id}{x}}f_s(t(\mathrm{}_0),h(\mathrm{}_0),u(\mathrm{}_0),N^1(\mathrm{}_0)).$$
(37)
Note that we do not know the ratio $`z_I/x`$, however (37) gives the form we should observe. More succinctly, for large $`N`$, this result predicts that for small $`t`$ and $`h`$

$$\mathrm{ln}Z_N(t,h)N^{1d}[\mathrm{ln}N]^p.$$
(38)
There is related work in number theory by Kanemitsu Kan  (cf. also Jap ). This paper studies moments of neighboring Farey fraction differences, which are similar to the “Farey tree” partition function F . At $`h=0`$, the latter has the same thermodynamics as the FFSC FK . However, Kan  uses a definition of the Farey fractions that, at each level, gives a subset of the Farey fractions employed here, and none of the moments considered corresponds to $`\beta =2`$ (the point of phase transition). It is interesting that, despite these differences, terms logarithmic in $`N`$ appear. More specifically, the sum of $`m`$th (integral) moments of the differences goes as

$$\frac{(\mathrm{ln}N)^{\delta _{2,m}}}{N^m}+O\left(\frac{(\mathrm{ln}N)^{h(m)}}{N^{m+g(m)}}\right),$$
(39)
for $`m2`$, with $`g(2)=1`$, $`g(3)=2`$ and $`g(m)=3`$ for $`m4`$, and $`h(m)=1`$ for $`2m4`$, $`h(m)=0`$ for $`m5`$. Now if all the Farey fractions were included (39) would apply to the Farey tree partition function with $`\beta =2m`$ (cf. FK ; F ) so that $`m2`$ would correspond to $`\beta 4`$. It would be interesting to determine whether (39) applies to the Farey tree partition function despite this difference, or to extend (39) to $`m=1`$ to see if it is consistent with (38).

## V 1-D KDP model with nonzero external field

In this section we consider the one-dimensional KDP (Potassium dihydrogen phosphate) model introduced by Nagle N . This model’s thermodynamics and energy level structure are similar to the Farey fraction spin chain, but it is easily solvable. Comparison of the two models thus sheds some light on the FFSC.

The KDP model exhibits first-order phase transitions only. The origin of the phase transition is infinite rather than long-range interactions.

The one-dimensional geometry of the model is illustrated in Fig. 3. It consists of $`N`$ cells, and each cell contains two dots. Each dot represents a proton in a hydrogen bond in the KDP molecule. Dots can be on the left or the right side of a cell. The energy of a neighboring pair of cells depends on the arrangement of dots at their common boundary. Only configurations with exactly two dots at each boundary (e.g. A, B and D in Fig. 3) are allowed, any other configuration (e.g. C in Fig. 3) has (positively) infinite energy and is therefore omitted. Of the allowed configurations, only two energies occur, $`0`$ (when there are two dots on the same side of a boundary, as in Fig. 3 D) or $`ϵ`$ (when the dots are on opposite sides, as in Fig. 3 A or B).

Let there be $`N`$ cells in a chain with periodic boundary conditions. Then there are two kinds of configurations with finite energy. In the first type of configuration, each cell has two dots on the same side. There are two such configurations and the total energy of each is $`0`$. In the second type of configuration, each cell has one dot on the left and one on the right. There are $`2^N`$ such configurations and the total energy of each is $`Nϵ`$. Thus, the partition function is simply

$$Z_N(\beta )=2+2^N\mathrm{exp}(\beta Nϵ).$$
(40)
It follows immediately that $`f=0`$ for $`\beta ϵ>\mathrm{ln}2`$ and $`f=ϵ\frac{\mathrm{ln}2}{\beta }`$ for $`\beta ϵ<\mathrm{ln}2`$. Thus the temperature of the (first-order) phase transition is $`T_c=ϵ/(\mathrm{ln}2)`$ and there is a latent heat with entropy change $`\mathrm{\Delta }s=\mathrm{ln}2`$. Clearly, the phase transition mechanism is a simple entropy-energy balance. At low temperatures, the ground state energy gives the minimal free energy, while in the high-temperature phase the extra entropy of the additional states gives a lower free energy.

Next, define the magnetization $`m`$ as the number of sides of cells with both dots on one side divided by the number of cells $`N`$. Then $`m=1`$ for $`\beta >\beta _c`$ and $`m=0`$ for $`\beta <\beta _c`$ (so that $`\mathrm{\Delta }m=1`$ at the phase transition), just as in the FFSC model.

Following the above definition of the magnetization, we introduce an external field $`h`$ by adding an energy $`\pm h/2`$ to each dot, according to whether it is on the right or left side of the cell. This gives the extra energy of an external field acting along the chain. Then the new partition function has the form

$$Z_N(\beta ,h)=2\mathrm{cosh}(\beta Nh)+2^N\mathrm{exp}(\beta Nϵ).$$
(41)
In the ordered phase $`Z\mathrm{exp}(\pm \beta Nh)`$. Thus, the free energy $`f=h`$, where the plus sign is for $`h>0`$ and the minus sign for $`h<0`$, exactly as in the FFSC. For the high-temperature phase $`Z\mathrm{exp}\left[N(\mathrm{ln}2\beta ϵ)\right]`$ and we get the same free energy as when $`h=0`$, $`f=ϵ\frac{\mathrm{ln}2}{\beta }`$. The phase boundary is given by $`h=\pm ϵt`$ (see Fig. 4), where $`t=\frac{\beta _c}{\beta }1`$ as before. Note the resemblance to the FFSC phase diagram (Fig. 2). Here, as $`\beta ϵ\mathrm{ln}2`$, $`h0`$ as it should, while for $`\beta 0`$ the field $`h\mathrm{ln}2/\beta `$. The entropy per site vanishes everywhere in the ordered phase, while for the high-temperature phase $`s=\mathrm{ln}2`$. Thus, this model has a non-zero latent heat and the phase transition is first-order everywhere. Note that the change in magnetization is $`\mathrm{\Delta }m=1`$ everywhere along the phase boundary between the ordered state and the high-temperature state.

Now for $`h=0`$, the FFSC has two ground states with all spins up or all spins down and energy independent of length $`N`$, just as in the KDP model. Then, in addition, the FFSC has $`2^N2`$ states with energies between $`\mathrm{ln}N`$ and $`Nc`$, for some constant $`c`$. On the other hand, the KDP model has just one energy ($`Nϵ`$) for the $`2^N`$ states corresponding to the $`2^N2`$ states of the Farey model. This might suggest that the states with energies close to $`\mathrm{ln}N`$ are responsible for the logarithmic factor in the Farey free energy, and thus shift the phase transition from first to second-order (for $`h=0`$). For $`h0`$ the energy of the $`\mathrm{ln}N`$ states is shifted by the field $`h`$ to order $`N`$, and the phase transition becomes first-order. However the mechanism of the FFSC phase transition may be more subtle. The “density of states” (number of configurations with a given energy) for the FFSC not well-behaved. In fact it is known rigorously that this quantity, summed over all chain lengths, has a limit distribution Pe .

Note that the free energy just derived is independent of $`h`$ in the high-temperature phase. Since this is not what we found for the FFSC, we consider another way to introduce an external field $`h`$ into the KDP model.

As before we have four different states for each cell. We index them with spin-one variables $`t_i`$ and $`s_i`$ ($`s_i,t_i\{0,+1,1\}`$ ) in each cell as in Fig. 5. Then the energy (for $`h=0`$) can be written

$$H_0=ϵ\underset{i=1}{\overset{N1}{}}t_i^2t_{i+1}^2$$
(42)
(assuming, in the sum, that the infinite energy contributions are omitted). The conditions $`s_i+t_i=\pm 1`$ and $`s_it_i=0`$ define the allowed states. We define the magnetization per site as

$$m=\frac{1}{N}\underset{i=1}{\overset{N}{}}(s_i+t_i).$$
(43)
Note that this definition gives a positive (negative) contribution if the upper dot in a given cell is on the right (left). (Note also that $`m^2=\frac{1}{N^2}_{i=1}^N(s_i+t_i)^2+\frac{1}{N^2}_{ij}(s_i+t_i)(s_j+t_j)=\frac{1}{N}_{j=2}^N(s_1+t_1)(s_j+t_j)+1/N`$.) Hence we can include an external field as follows

$$H=H_0h\underset{i=1}{}(s_i+t_i)=H_0hNm.$$
(44)
Thus

$$Z(\beta ,h)=e^{\beta Nh}+e^{\beta Nh}+e^{\beta ϵN}[2\mathrm{cosh}(\beta h)]^N,$$
(45)
and the free energy in high-temperature phase becomes

$$f(\beta ,h)=ϵ\frac{\mathrm{ln}(2\mathrm{cosh}(\beta h))}{\beta }$$
(46)
or for small $`h`$

$$ftϵ\frac{\mathrm{ln}2}{2ϵ(t+1)}h^2,$$
(47)
with $`t=\frac{\beta _c}{\beta }1`$ as above. The phase boundary is given by

$$\beta h=\mathrm{ln}\left(2\mathrm{cosh}(\beta h)\right)\beta ϵ.$$
(48)
For $`\beta h1`$ and $`h>0`$, using $`\beta _c=\frac{\mathrm{ln}2}{ϵ}`$, this gives

$$h=ϵt+\frac{ϵ\mathrm{ln}2}{2}t^2+O(t^3),$$
(49)
The phase diagram near the critical point is very close to the previous one (see Fig. 4). The magnetization in the ordered phase is again independent of temperature, i.e. $`m=\pm 1`$. In the high-temperature phase we have $`m=\mathrm{tanh}(\beta h)`$. Thus the magnetization change across the phase boundary close to the critical point is $`\mathrm{\Delta }m=1t\mathrm{ln}2`$. The transition is again first-order, with the entropy change $`\mathrm{\Delta }s=\mathrm{ln}2(1\frac{\mathrm{ln}2}{2}t^2)`$. Results for $`h<0`$ follow immediately by symmetry.

## VI Summary and comments

In this paper, we have extended the definition of the Farey fraction spin chain to include an external field $`h`$. From rigorous and more heuristic arguments, we have determined the phase diagram and phase transition behavior of the extended model. Our results are fully consistent with scaling theory (for the case when a “marginal” field is present) despite the unusual nature of the transition for $`h=0`$. In particular, we find for the renormalization group eigenvalues $`y_h=y_t=d`$, and for the sub-leading eigenvalues $`z_t=x`$ and $`z_h=0`$. We also examine a completely solvable model with very similar thermodynamics, but for which all phase transitions are first-order.

## VII Acknowlegements

We are grateful to J. L. Cardy and M. E. Fisher for useful suggestions. We also thank an anonymous referee for alerting us to an incompleteness in one of our proofs. This work was supported in part by the National Science Foundation Grant No. DMR-0203589. *

## Appendix A Bounds for $`\frac{Z_{N+1}}{Z_N}`$

First we introduce some notation (following FK ).

We use $`r_N^{(n)}:=\frac{n_N^{(n)}}{d_N^{(n)}}`$ for the fractions (called Farey fractions), where $`n`$ is the order of the Farey fraction in level $`N`$. Level $`N=0`$ consists of the two fractions $`\{\frac{0}{1},\frac{1}{1}\}`$. Succeeding levels are generated by keeping all the fractions from level $`N`$ in level $`N+1`$, and including new fractions. The new fractions at level $`N+1`$ are defined via $`d_{N+1}^{(2n)}:=d_N^{(n)}+d_N^{(n+1)}`$ and $`n_{N+1}^{(2n)}:=n_N^{(n)}+n_N^{(n+1)}`$, so that
$`N=0\{\frac{0}{1},\frac{1}{1}\}`$
$`N=1\{\frac{0}{1},\frac{1}{2},\frac{1}{1}\}`$
$`N=2\{\frac{0}{1},\frac{1}{3},\frac{1}{2},\frac{2}{3},\frac{1}{1}\}`$, etc.
Note that $`n=1,\mathrm{},2^N+1`$. When the Farey fractions are defined using matrices (spin states) A and B, the level $`N+1`$ is the number of matrices in the chains starting with matrix $`A`$ and hence the length of the spin chain K-O .

Using this notation we can write the partition function (3) restricted to chains starting with $`A`$

$$Z_N^A(\beta ,h)=\underset{n=1}{\overset{2^N}{}}\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)}}{(d_N^{(n)}+n_N^{(n+1)})^\beta },\beta .$$
(50)
Note that the partition function (3) is the sum of $`Z_N^A(\beta ,h)`$ and $`Z_N^B(\beta ,h)`$, where the $`Z_N^B(\beta ,h)`$ is the partition function for chains starting with the matrix $`B`$. First we find bounds for $`Z_N^A(\beta ,h)`$ and then prove a lemma which lets us apply the bounds for $`Z_N^A(\beta ,h)`$ to $`Z_N^B(\beta ,h)`$ also.

Now, when we go from level $`N`$ to level $`N+1`$ we double the number of the terms in the partition function. Note that for chains starting with the matrix $`A`$ one half of the terms come from matrix products of the form $`AM_{N1}A`$ and the others from products $`AM_{N1}B`$. It is easy to check that the corresponding traces for given $`n\{1,\mathrm{},2^N\}`$ are $`d_{N+1}^{(2n1)}+n_{N+1}^{(2n)}`$ and $`d_{N+1}^{(2n)}+n_{N+1}^{(2n+1)}`$, respectively. These traces are multiplied by an $`h`$ dependent factor $`e^{\beta h\left(2_{i=1}^{N+1}\sigma _iN1\right)}`$ which is simply $`e^{\beta h}`$ raised to the power $`(\mathrm{\#}A\mathrm{\#}B)`$, the number of matices $`A`$ minus the number of matices $`B`$ in the particular chain. For the terms from products of the form $`AM_{N1}A`$, it follows on using the definition of the Farey fractions that

$$\frac{e^{\beta h\left(2_{i=1}^{N+1}\sigma _iN1\right)}}{(d_{N+1}^{(2n1)}+n_{N+1}^{(2n)})^\beta }=\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)+\beta h}}{(d_N^{(n)}+n_N^{(n)}+n_N^{(n+1)})^\beta }\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)}}{(d_N^{(n)}+n_N^{(n+1)})^\beta }e^{\beta |h|}$$

and, similarly, for $`AM_{N1}B`$

$$\frac{e^{\beta h\left(2_{i=1}^{N+1}\sigma _iN1\right)}}{(d_{N+1}^{(2n)}+n_{N+1}^{(2n+1)})^\beta }=\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)\beta h}}{(d_N^{(n)}+d_N^{(n+1)}+n_N^{(n+1)})^\beta }\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)}}{(d_N^{(n)}+n_N^{(n+1)})^\beta }e^{\beta |h|}.$$

For the lower bound we just need the $`AM_{N1}A`$ terms

$$\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)+\beta h}}{(d_N^{(n)}+n_N^{(n)}+n_N^{(n+1)})^\beta }\frac{e^{\beta h\left(2_{i=1}^N\sigma _iN\right)}}{(d_N^{(n)}+n_N^{(n+1)})^\beta }\frac{e^{\beta |h|}}{2^\beta },$$

where we used the fact that $`n_N^{(n)}d_N^{(n)}`$. Thus we get for $`Z_{N+1}^A(\beta ,h)=Z_{N+1}^{AM_{N1}A}(\beta ,h)+Z_{N+1}^{AM_{N1}B}(\beta ,h)`$

$$2^\beta e^{\beta |h|}Z_N^A(\beta ,h)Z_{N+1}^A(\beta ,h)2e^{\beta |h|}Z_N^A(\beta ,h)$$

for any $`\beta 0`$ and $`h.`$

Finally, we prove a lemma which allows us to bound $`Z_N^B(\beta ,h)`$. Consider a $`(2\times 2)`$ matrix $`M=\left(\genfrac{}{}{0pt}{}{m_1}{m_3}\genfrac{}{}{0pt}{}{m_2}{m_4}\right)`$ and define the operator $``$ via $`\stackrel{~}{M}:=\left(\genfrac{}{}{0pt}{}{m_4}{m_2}\genfrac{}{}{0pt}{}{m_3}{m_1}\right)`$. Then we have the following result.

###### Lemma A.1

Let $`M=AZ_1Z_2\mathrm{}Z_N`$, where $`Z_i\{A,B\}`$, with $`A=\left(\genfrac{}{}{0pt}{}{1}{1}\genfrac{}{}{0pt}{}{0}{1}\right)`$ and $`B=\left(\genfrac{}{}{0pt}{}{1}{0}\genfrac{}{}{0pt}{}{1}{1}\right)`$. Then $`\stackrel{~}{M}=B\stackrel{~}{Z}_1\stackrel{~}{Z}_2\mathrm{}\stackrel{~}{Z}_N`$, i.e. the $``$ operator exchanges $`A`$ and $`B`$.

Proof. We will use mathematical induction. It is easy to see that $`A=\stackrel{~}{B}`$ and $`B=\stackrel{~}{A}`$. From matrix multiplication follows $`B\stackrel{~}{M}=\left(\genfrac{}{}{0pt}{}{m_2+m_4}{m_2}\genfrac{}{}{0pt}{}{m_1+m_3}{m_1}\right)`$ and $`AM=\left(\genfrac{}{}{0pt}{}{m_1}{m_1+m_3}\genfrac{}{}{0pt}{}{m_2}{m_2+m_4}\right)`$.

Clearly the $``$ operation is a 1-to-1 map of the set of all chains $`AM_N`$ onto $`BM_N`$. Furthermore, the magnetic field term in the energy of each chain changes sign under this operation, so that the bounds just obtained for $`Z_N^A(\beta ,h)`$ may be applied to $`Z_N^B(\beta ,h)`$. Therefore

$$2^\beta e^{\beta |h|}\frac{Z_{N+1}}{Z_N}2e^{\beta |h|}.$$

Note that the proof is easily adapted to the KSC model.