format: text

text:
# Generalized Fermi-Dirac Functions and Derivatives: Properties and Evaluation

## 1 Introduction

Relativistic quantum statistics is important for studying thermodynamic properties of electrons at the elevated degeneracy and relatively high temperature as found in the interior of massive stars. Chandrasekhar Chandra  expresses the thermodynamic quantities as integrals of hyperbolic functions \[see Eq. (170) - (172) of Chandra \]. In a more general way Cox & Giuli Cox , write them in term of the generalized Fermi-Dirac (FD) functions with different orders \[see Eq. (24.98) - (24.100) of Cox \]. They define the generalized FD functions as

$$F_k(\eta ,\beta )=_0^{\mathrm{}}\frac{x^k\sqrt{1+(\beta x/2)}}{\mathrm{}^{x\eta }+1}x,$$
(1)
where $`\eta \mu /k_\mathrm{B}T`$ and $`\beta k_\mathrm{B}T/m_ec^2`$ are the degeneracy parameter and a dimensionless temperature, respectively. Here, $`m_e`$ is electron mass, $`T`$ absolute temperature, $`\mu `$ the chemical potential, $`c`$ speed of light in vacuum and $`k_\mathrm{B}`$ the Boltzmann constant. The non-relativistic form of Eq. (1) (corresponding to $`\beta =0`$) is

$$F_k(\eta )=_0^{\mathrm{}}\frac{x^k}{\mathrm{}^{x\eta }+1}x$$
(2)
is usually called the Fermi-Dirac function, and is used to study thermodynamic properties of the degenerate non-relativistic electron gas. For the usual physical applications, both for the non-relativistic and relativistic functions, the index $`k`$ is an integer or half integer with $`k1`$.

No exact analytic formulae are available to evaluate the generalized FD functions in full range of degeneracy and temperature. Some approximate expressions exist for some extreme cases, such as when degeneracy is very low, or extremely high Cox ; Mira ; Pichon ; Blud . There are also asymptotic expressions for the thermodynamic quantities alone Blin . On the other hand, Eggleton et al. Egg , give a very useful global fitting formula for the generalized FD functions \[see also Pichon \], even though its accuracy is relatively modest.

Numerical evaluation of the generalized FD functions is not an easy task. Simple tabular interpolations Cox  can cover only a small range of degeneracy and temperature, and the accuracy is not very good. Because the integrand is proportional to $`x^k`$ when $`x0`$, and decays as $`\mathrm{}^x`$ when $`x\mathrm{}`$, few simple numeric integration methods can simultaneously cover a sufficiently large range of degeneracy and temperature with acceptable accuracy as well as with an acceptable calculation effort. Sagar  uses a modified Gauss quadrature with weight function $`w(x)=x^k/(\mathrm{}^x+1)`$, and Gau  uses weight function $`w(x)=x^k/\mathrm{}^x`$ to integrate the generalized FD functions directly. Both methods are relatively fast, but are only good for $`\eta 5`$. In Pichon , a Gauss-Laguerre quadrature with adjusted decay factor is used to cover $`1\eta 11`$, and approximation formulae beyond that region. In Mira ; NaMo , different variations of the trapezoidal rule are being used to do the integration, but both involve a heavy trade-off between accuracy and computing time.

Mohankumar & Natarajan MoNa  use another approximation method by subtracting the sum of residues of finite number of poles within a pre-selected contour from results of trapezoidal integration. This method is fast when $`\eta `$ is not too big, and the results can reach 14-digit accuracy in double precision for $`\beta <100`$, and when $`\eta `$ is not too small, typically from $`\eta >25`$ to $`\eta =2000`$. It is good to about 10 digits when $`\eta `$ is as large as 10000. However, the number of poles to be included increases rapidly when $`\eta `$ gets too big. Extending this method to calculate different derivatives of the generalized FD functions with the same order of accuracy is not an easy task.

Aparicio Apa  proposed another method which consists of splitting the generalized-FD-function integral into four regions. The break points are chosen such that errors from all four regions are balanced, and the location of the break points is optimized. Based on the property of the integrand, the Gauss-Legendre quadratures used are from $`x=0`$ to $`x=s_3`$, and the Gauss-Laguerre quadrature from $`x=s_3`$ to $`x\mathrm{}`$ \[see Appendix (B)\]. The method is fast (proportional to the number of points used in Gauss quadratures) and accurate. Extending it to include derivatives of the generalized FD functions is straightforward. Our numerical method to evaluate derivatives of the generalized FD functions follows this scheme closely. We discuss it in more detail in Sect. 3.

Besides the generalized FD functions themselves, their derivatives are also important for computing the thermodynamic quantities of an electron gas. Miralles & Van Riper Mira  give the results of the first order of derivatives of the generalized FD functions with respect to both $`\eta `$ and $`\beta `$. They use approximation formulae when $`\eta `$ is either very large or very small, and interpolate through tables when $`\eta `$ is in between. The disadvantage is that it is only accurate to about 6 digits, and discontinuity exists when two different methods join. Also, their results are limited to $`k=\frac{1}{2},\frac{3}{2}`$ and $`\frac{5}{2}`$ only. On the other hand, solar equation of state studies require up to third order derivatives on the generalized FD functions CD92 . In particular, the precise data from helioseismology have demonstrated that the relativistic electron gas is necessary for solar modeling EK98 . For solar applications, high-precision computations of the relativistic thermodynamic quantities are mandatory Gong .

Because of its complexity, the generalized FD functions have not been fully understood. In Sect. 2 we first present some newly found recursion relations, and then discuss some of their properties. The numerical method to calculate derivatives of the generalized FD functions up to the third order is introduced in Sect. 3, based on the scheme of Apa . At the end we check the efficiency and validity of some popularly used analytic approximation formulae in several extreme cases in Sect. 4.

## 2 Recursion Relations

Cox & Giuli Cox  show that the generalized FD functions obey the following recursion relations

$$\beta \frac{F_k(\eta ,\beta )}{\beta }=\frac{F_{k+1}(\eta ,\beta )}{\eta }(k+1)F_k(\eta ,\beta )$$
(3)
which link the FD functions $`F_k(\eta ,\beta )`$ and $`F_{k+1}(\eta ,\beta )`$ via both their $`\eta `$ and $`\beta `$ derivatives. Beyond that, we have found that the generalized FD functions also obey another recursion relation

$$F_{k+1}(\eta ,\beta )=4\frac{F_k(\eta ,\beta )}{\beta }+2\beta \frac{F_{k+1}(\eta ,\beta )}{\beta },$$
(4)
which is between the generalized FD functions and their $`\beta `$ derivatives.

From Eq. (3) and Eq. (4) we can further derive two useful recursion relations, in which

$$\frac{F_{k+1}(\eta ,\beta )}{\eta }=\left(k+\frac{3}{2}\right)F_k(\eta ,\beta )2\frac{F_{k1}(\eta ,\beta )}{\beta }$$
(5)
is a relation among $`F_{k1}(\eta ,\beta )`$, $`F_k(\eta ,\beta )`$ and $`F_{k+1}(\eta ,\beta )`$, with both $`\eta `$ and $`\beta `$ derivatives, and

$$kF_{k1}(\eta ,\beta )+\left(\frac{1}{2}k+\frac{3}{4}\right)\beta F_k(\eta ,\beta )=\frac{F_k(\eta ,\beta )}{\eta }+\frac{1}{2}\beta \frac{F_{k+1}(\eta ,\beta )}{\eta }$$
(6)
is between the generalized FD functions and their $`\eta `$ derivatives only.

These recursion relations can be used to check the accuracy of numerical evaluations as well as to simplify some theoretical expressions. However, one should keep in mind that when $`\eta `$ is very small the two terms on the right-hand side of either Eq. (3) and Eq. (5) are so close to each other that almost perfect cancellation happens. For this reason, the use of recursion relations to evaluate the derivatives of the generalized FD functions is not always possible.

One interesting property of the recursion relation in Eq. (3) is that, for the non-relativistic limit case ($`\beta =0`$), it can be written as:

$$\frac{F_{k+1}(\eta )}{\eta }=(k+1)F_k(\eta )$$
(7)
which can be further written as:

$$_0^{\mathrm{}}\frac{x^k}{\mathrm{}^{x\eta }+1}x=_0^{\mathrm{}}\frac{x^k}{\mathrm{}^{x\eta }+1}\frac{x}{(\mathrm{}^{x\eta }+1)(k+1)}x$$
(8)
Both Eqs. (7) and (8) are satisfied for arbitrary $`\eta `$. This is very remarkable, because although the function $`x/[(\mathrm{}^{x\eta }+1)(k+1)]`$ is not equal to $`1`$, when multiplied with the Fermi-Dirac kernel $`x^k/(\mathrm{}^{x\eta }+1)`$, it leads to exactly the same value of the definitive integrals as if it were equal to $`1`$.

Another interesting behavior we observed from the generalized FD functions is the dependence on $`\eta `$ of the following ratio (see Fig. 1).

$$\frac{\{[F_k(\eta ,\beta )/\eta ]+\beta [F_{k+1}(\eta ,\beta )/\eta ]\}^2}{[F_k(\eta ,\beta )+\beta F_{k+1}(\eta ,\beta )]\{[^2F_k(\eta ,\beta )/\eta ^2]+\beta [^2F_{k+1}(\eta ,\beta )/\eta ^2]\}}$$
(9)
The expression appears in the derivatives of number density or internal energy of the relativistic electron gas (see GD  for more detail).

For small degeneracy ($`\eta \mathrm{}`$) the ratio (9) is very close to 1 for all value of index $`k`$ ($`k>1`$) and all $`\beta 0`$. However, when the electron gas is extremely degenerate ($`\eta \mathrm{}`$), the ratio (9) approaches different limits for the non-relativistic and relativistic cases. For the non-relativistic case ($`\beta =0`$), the limits of (9) are $`\frac{1}{1},\frac{3}{1},\frac{5}{3},\frac{7}{5},\mathrm{}`$ when $`k=\frac{1}{2},\frac{1}{2},\frac{3}{2},\frac{5}{2},\mathrm{}`$, respectively. The general expression is $`\frac{k+1}{k}`$ for $`k>1`$ and $`k0`$. For the relativistic case ($`\beta 0`$) the limits of (9) are $`\frac{2}{1},\frac{3}{2},\frac{4}{3},\frac{5}{4},\mathrm{}`$ when $`k=\frac{1}{2},\frac{1}{2},\frac{3}{2},\frac{5}{2},\mathrm{}`$, respectively. The general expression is $`\frac{k+5/2}{k+3/2}`$ ($`k>1`$). This result is true for $`\eta \beta 1`$. So for a given small $`\beta `$ ($`\beta 0`$), that is, in the nearly non-relativistic case, for sufficiently large $`\eta `$, the relativistic relation is satisfied, as shown in Fig. (1). This property links the behavior of the generalized FD functions of both extreme cases (very small degeneracy or extremely great degeneracy). In Appendix (D) we will give a proof of this property using asymptotic expansions.

One thing we want to point out is that, because the ratio of Eq. (9) is very close to unity when $`\eta `$ is small, one has to be cautious when evaluating the difference between the nominator and denominator for the case of a slightly degenerate electron gas in order to avoid numerical noise due to near cancellation.

## 3 Numerical Method and Test Run Results

We choose the scheme proposed by Apa  to calculate derivatives of the generalized FD functions, because this method is fast, accurate and transparent. The expressions for the derivatives are listed in Appendix (A). In the method adopted, the integration is further split into four domains. Gauss-Legendre quadratures are applied to the first three pieces, starting from $`x=0`$ to a preselected point $`s_3`$, as defined in Appendix (B). In the last part, where $`x`$ is from $`s_3`$ to $`\mathrm{}`$ the $`\mathrm{}^x`$ term dominates, and hence Gauss-Laguerre quadrature is used. The positions of the break points are optimized such that the combined errors from all the pieces contributing to the generalized FD functions are minimized and equally distributed among the pieces. Fortunately, the choice of break points turns out rather uncritical for each individual derivative \[see also Apa \].

As shown in Appendix (A), the integrand of each derivative can be represented as a product of the generalized FD kernel $`\left[x^k\sqrt{1+(\beta x/2)}\right]/(\mathrm{}^{x\eta }+1)`$ with a tuning function. When $`\eta `$ and $`\beta `$ are not very large, this tuning function varies mildly \[see Fig. (2)\]. Hence the same break points chosen for the generalized FD functions themselves can also be applied to their derivatives, and good accuracy is achieved. However, when $`\eta `$ and $`\beta `$ are getting large, new break points must be chosen. The $`\beta `$ derivatives are not very sensitive to the choice of break points. Therefore, we provide new break points only for $`\eta `$ derivatives. Formulae to calculate break points can be found in Appendix (B).

One concern in the calculation of the generalized FD functions and their derivatives is how one determines the accuracy of the results. For the generalized FD functions we compare the results of Apa  (the method adopted in the present paper) with those of MoNa  evaluated in double precision. We find that they agree with each other to 14-digit accuracy for $`\eta `$ up to 1000, and to 10-digit accuracy for $`\eta `$ up to 10000. Additionally, both of them agree with Mira  to 8-to-10-digit accuracy for $`\eta `$ up to 1000.

For the derivatives of the generalized FD functions, we first test them using the different recursion relations available. Except for the aforementioned cancellations when $`\eta `$ is small, the recursion relations are satisfied with better than 10-digit accuracy, for all derivatives with $`\beta `$ being between $`10^6`$ and $`10^4`$ and $`\eta `$ between -100 and 10000. The only exception is the third-order derivative with respect to $`\eta `$. Here, the same accuracy can only be achieved for $`\eta `$ less than 500. All of the above results refer to double precision (128 bits) (which is single precision on a Cray). Of course, better results can be achieved with quadruple precision. In practice, if $`k\frac{1}{2}`$, the recursion relations can be used to convert the third-order $`\eta `$ derivative into a lower-order one when $`\eta `$ is big, and in this way, very accurate results can be achieved.

As a second test, we compare our directly evaluated derivatives with those obtained from numerical differentiation. For example, we compare $`F_{1/2}(\eta ,\beta )/\eta `$ from our method with $`\left[F_{1/2}(\eta +\mathrm{\Delta }\eta ,\beta )F_{1/2}(\eta \mathrm{\Delta }\eta ,\beta )\right]/2\mathrm{\Delta }\eta `$, with the difference $`\mathrm{\Delta }\eta `$ chosen to be $`10^2`$, $`10^4`$, $`10^6`$, $`10^8`$, respectively, and the calculation is done in double precision. The purpose of the test is (i) to test the accuracy of our method, and (ii) to find the dependence of $`\mathrm{\Delta }\eta `$ as a function of accuracy.

We conclude that both methods agree well. In practice, a step of $`\mathrm{\Delta }\eta `$ or $`\mathrm{\Delta }\beta `$ between $`10^2`$ and $`10^4`$ is an optimal choice to compute the derivatives of the generalized FD functions by numerical differentiation. For example, the numerical evaluation of the second-order derivative with respect to $`\eta `$ from the first-order derivative yields an accuracy of at least 8 digits in double precision.

However, we would like to stress that our direct evaluations are superior to the numerical differentiation used here only for testing purposes. The accuracy of numerical differentiation can drop to half of the machine precision \[i.e. to 8 to 9 digits, compared to the 16 digit accuracy of double precision (non-Cray machines) or single precision (Cray machines)\] when the aforementioned 3-point numerical differentiation is applied once. Applied twice, the resulting accuracy can drop to 3 to 5 digits for $`\eta <2000`$, and only 2 digits when $`\eta 5000`$. Applied three times, the resulting accuracy would only be 2 digits when $`\eta 50`$, and for higher $`\eta `$ values, the results would be completely useless. This demonstrates the well-known fact that numerical differentiation is a dangerous procedure if it is used beyond first-order derivatives. In addition, the step size to be used has to be chosen judiciously. For high-precision applications such as found in helioseismology, which requires precise second- and third-order thermodynamic quantities, the direct approach method presented in this paper is by far superior to more-than-one-time numerical differentiations.

In this paper we have always used 200-point Gauss quadratures. Although fewer-points quadratures are also good for the calculations of the generalized FD functions Apa , the 200-point Gauss quadratures used here guarantee the accuracy of the high-order derivatives. If computing speed is important, and third-order $`\eta `$-derivatives are not needed, then the 40-to-80-point Gauss quadratures in Apa  are a good choice.

To use the program accompanying this paper, the user has first to compile and link the three Fortran files. Then, upon choosing the index of the requested generalized FD function $`k`$, the label $`\mathrm{IB}`$ denoting either the function itself or a specified derivative (see Table 1), the range of degeneracy $`\eta `$ and dimensionless temperature $`\beta `$, the program then calculates the desired generalized FD function or derivative. A sample input is given in the file input, and the corresponding results are contained in the file testResult.

## 4 Validity of Analytic Approximations

Although the present numerical method to evaluate generalized FD functions is accurate within a very large range of input parameters $`\eta `$ and $`\beta `$, analytic approximations still have advantages for theoretical studies and the understanding of the underlying physics. However, since analytic approximations are usually only valid for restricted parameter regions, we must know them for each of them. In the following, we limit ourselves to the properties of the $`k=\frac{1}{2}`$ functions only if its behavior is representative for the other cases. Otherwise, we discuss the cases of index $`k=\frac{1}{2}`$ and $`k=\frac{3}{2}`$.

First we check the accuracy of the non-relativistic (NR) FD function as an approximation to the generalized FD functions in the weakly non-relativistic case ($`\beta 1`$)

$$F_k(\eta ,\beta )F_k(\eta ).$$
(10)
From Fig. 3, we can see that when both $`\eta `$ and $`\beta `$ are small the approximation is good to 6 digits, and the accuracy drops when either $`\eta `$ or $`\beta `$ is getting big. Considerably better agreement (twice as many accurate digits, as shown in Fig. 4) can be reached by including one more term is Blud :

$$F_k(\eta ,\beta )F_k(\eta )+\frac{1}{4}\beta F_{k+1}(\eta )$$
(11)
For the ultra-relativistic case ($`\beta 1`$) Blud , results of the following approximation is shown in Fig. 5:

$$F_k(\eta ,\beta )\sqrt{\frac{\beta }{2}}\left[F_{k+\frac{1}{2}}(\eta )+\frac{1}{\beta }F_{k\frac{1}{2}}(\eta )\right]$$
(12)
The figure shows that more than 6-digit accuracy can be obtained with this formula if $`\beta >0`$, and $`\eta `$ is not too small, typically $`\eta 0`$, and more than 10-digit accuracy when $`\beta >2`$.

When $`\eta >0`$ we can see that the asymptotic expansions by Pichon Pichon  \[see Appendix C.1, Eq. (C1-C11); also Cox & Giuli Cox \] are very good for $`F_{\frac{1}{2}}(\eta ,\beta )`$ (Fig. 6), with better than 12-digit accuracy when $`\eta >120`$. However, the results are only moderate for $`F_{\frac{3}{2}}(\eta ,\beta )`$ (Fig. 7) and $`F_{\frac{5}{2}}(\eta ,\beta )`$, in which case the accuracy is only better than 6 digits when $`\eta >5000`$, and much less accurate when $`\eta `$ is smaller. For convenience, We list the formulae for large and small $`\eta `$ in Appendix (C).

When $`\eta `$ is small, the expression of the generalized FD functions in terms of modified Bessel functions \[see Appendix C.2, Eq. (C12-C17); Cox & Giuli Cox \] achieve very high accuracy. We have used Mathematica 4.0 (for Solaris) to evaluate the modified Bessel functions in double precision. Since the leading terms in Eq. (43) - (45) are comparable when $`\beta `$ is small, cancellation effect would appear. A power series expansion expression Mira  is therefore more appropriate. We have combined the methods and show their results in Fig. 8.

We have also checked the quality of the fitting formulae by Egg  \[see Appendix C.3, Eq. (C21-C26)\] (Fig. 9). Specifically, we have adopted the expression given by Pichon  (his Eq. 5a - 5c). These fitting formulae cover a much larger range of $`\eta `$ and $`\beta `$ than those asymptotic formulae above (but with lower accuracy), and they are adequate if accuracy is not critical, and temperature not too high. For instance, with these fitting formulae we can obtain 4-digit accuracy when $`\eta <1000`$ and $`\beta <10^4`$, and 2-digit accuracy when $`\eta <10000`$ and $`\beta <3\times 10^5`$.

## 5 Conclusions

In this paper we have investigated properties of the generalized FD functions in some detail, and we have found some useful new recursion relations. We propose an accurate and relatively fast numerical method to evaluate the derivatives of the generalized FD functions, following the algorithm introduced by Apa . Limitations of direct numerical differentiation are also discussed. Finally, we have obtained the range of validity of some popular analytic approximations for the generalized FD functions.

ZG would like to thank A. Natarajan and N. Mohabkumar for providing us with their code to evaluate generalized FD functions, and A. Nayfonov for useful discussions. ZG, LZ and WD are supported in part by the NSF grant AST-9618549 and AST-9987391.

## Appendix A APPENDIX: Derivatives of the Generalized FD Functions

Derivatives of the generalized FD functions up to third order with respect to $`\eta `$ and $`\beta `$.

$$\frac{F_k(\eta ,\beta )}{\eta }=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{1}{1+\mathrm{}^{\eta x}}x$$
(13)
$$\frac{F_k(\eta ,\beta )}{\beta }=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{x}{4+2\beta x}x$$
(14)
$$\frac{^2F_k(\eta ,\beta )}{\eta ^2}=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{1\mathrm{}^{\eta x}}{(1+\mathrm{}^{\eta x})^2}x$$
(15)
$$\frac{^2F_k(\eta ,\beta )}{\eta \beta }=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{x}{4+2\beta x}\frac{1}{1+\mathrm{}^{\eta x}}x$$
(16)
$$\frac{^2F_k(\eta ,\beta )}{\beta ^2}=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{x^2}{(4+2\beta x)^2}x$$
(17)
$$\frac{^3F_k(\eta ,\beta )}{\eta ^3}=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{(1\mathrm{}^{\eta x})^22\mathrm{}^{\eta x}}{(1+\mathrm{}^{\eta x})^3}x$$
(18)
$$\frac{^3F_k(\eta ,\beta )}{\eta ^2\beta }=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{1\mathrm{}^{\eta x}}{(1+\mathrm{}^{\eta x})^2}\frac{x}{4+2\beta x}x$$
(19)
$$\frac{^3F_k(\eta ,\beta )}{\eta \beta ^2}=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{x^2}{(4+2\beta x)^2}\frac{1}{1+\mathrm{}^{\eta x}}x$$
(20)
$$\frac{^3F_k(\eta ,\beta )}{\beta ^3}=_0^{\mathrm{}}\frac{x^k\sqrt{1+\frac{\beta x}{2}}}{\mathrm{}^{x\eta }+1}\frac{3x^3}{(4+2\beta x)^3}x$$
(21)
## Appendix B APPENDIX: Break Points

Three break points $`S_1,S_2,S_3`$ are \[as shown in Apa \]:

$`S_1`$ $`=`$ $`X_aX_b`$ (22)
$`S_2`$ $`=`$ $`X_a`$ (23)
$`S_3`$ $`=`$ $`X_a+X_c`$ (24)
where

$$X_a=\frac{a_1+b_1\xi +c_1\xi ^2}{1+c_1\xi }$$
(25)
$$X_b=\frac{a_2+b_2\xi +c_2d_2\xi ^2}{1+e_2\xi +c_2\xi ^2}$$
(26)
$$X_c=\frac{a_3+b_3\xi +c_3d_3\xi ^2}{1+e_3\xi +c_3\xi ^2}$$
(27)
$$\xi (\eta )=\sigma ^1\mathrm{ln}[1+\mathrm{}^{\sigma (\eta D)}]$$
(28)
where the corresponding parameters are:

| Parameter | $`F_k(\eta ,\beta )`$ | $`\frac{F_k(\eta ,\beta )}{\eta }`$ | $`\frac{^2F_k(\eta ,\beta )}{\eta ^2}`$ | $`\frac{^3F_k(\eta ,\beta )}{\eta ^3}`$ |
| --- | --- | --- | --- | --- |
| $`D`$ | 3.3609E0 | 4.99551E0 | 3.93830E0 | 4.17444E0 |
| $`\sigma `$ | 9.1186E-2 | 9.11856E-2 | 9.11856E-2 | 9.11856E-2 |
| $`a_1`$ | 6.7774E0 | 6.77740E0 | 6.77740E0 | 6.77740E0 |
| $`b_1`$ | 1.1418E0 | 1.14180E0 | 1.14180E0 | 1.14180E0 |
| $`c_1`$ | 2.9826E0 | 2.98255E0 | 2.98255E0 | 2.98255E0 |
| $`a_2`$ | 3.7601E0 | 3.76010E0 | 3.76010E0 | 3.76010E0 |
| $`b_2`$ | 9.3719E-2 | 9.37188E-2 | 9.37188E-2 | 9.37188E-2 |
| $`c_2`$ | 2.1064E-2 | 2.10635E-2 | 2.10635E-2 | 2.10635E-2 |
| $`d_2`$ | 3.1084E+1 | 3.95015E1 | 3.14499E1 | 3.05412E1 |
| $`e_2`$ | 1.0056E0 | 1.00557E0 | 1.00557E0 | 1.00557E0 |
| $`a_3`$ | 7.5669E0 | 7.56690E0 | 7.56690E0 | 7.56690E0 |
| $`b_3`$ | 1.1695E0 | 1.16953E0 | 1.16953E0 | 1.16953E0 |
| $`c_3`$ | 7.5416E-1 | 7.54162E0 | 7.54162E0 | 7.54162E0 |
| $`d_3`$ | 6.6559E0 | 7.64734E0 | 6.86346E0 | 7.88030E0 |
| $`e_3`$ | -1.2819E0 | -1.28190E-1 | -1.28190E-1 | -1.28190E-1 |

## Appendix C APPENDIX: Asymptotic Expressions for Great and Small Degeneracy

### C.1 Great degeneracy ($`\eta >0`$) Cox ; Pichon

$`F_{\frac{1}{2}}(\eta ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\beta ^3}}}\left\{f_{\frac{1}{2}}(y)+(1+\eta \beta )[C_1+C_2+C_3(4y^2+7)]\right\}`$ (29)
$`F_{\frac{3}{2}}(\eta ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\beta ^5}}}\{f_{\frac{3}{2}}(y)+C_1(3+2\eta \beta )C_2`$ (30)
$`C_3[(4\eta \beta +6)\eta \beta +3]\}`$
$`F_{\frac{5}{2}}(\eta ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\beta ^7}}}(f_{\frac{5}{2}}(y)+C_1(5+\eta \beta )+C_2\{[(2\eta \beta +10)\eta \beta +15]\eta \beta +5\}`$ (31)
$`+C_3(3+5\eta \beta ))`$
where

$$1+y^2=(1+\eta \beta )^2$$
(32)
$$C_1=\frac{\pi ^2}{6}\frac{\beta ^2}{y},C_2=\frac{7}{20}C_1\left(\frac{\pi \beta }{y^2}\right)^2,C_3=\frac{31}{168}C_1\left(\frac{\pi \beta }{y^2}\right)^4$$
(33)
For $`\eta 0.05`$:

$$f_{\frac{1}{2}}=\frac{1}{2}\left[y\sqrt{1+y^2}\mathrm{arcsinh}(y)\right]$$
(34)
$$f_{\frac{3}{2}}=\frac{1}{3}y^3f_{\frac{1}{2}}(y)$$
(35)
$$f_{\frac{5}{2}}=\frac{5}{8}y\left(1+\frac{2}{5}y^2\right)\sqrt{1+y^2}\frac{2}{3}y^3\frac{5}{8}\mathrm{arcsinh}(y)$$
(36)
For $`\eta <0.05`$:

$`f_{\frac{1}{2}}`$ $`=`$ $`(\{[(\{[({\displaystyle \frac{58773}{1114112}}y^2+{\displaystyle \frac{77}{5120}})y^2{\displaystyle \frac{63}{3328}}]y^2+{\displaystyle \frac{35}{1408}}\}y^2{\displaystyle \frac{5}{144}})y^2`$ (37)
$`+{\displaystyle \frac{5}{36}}]y^2{\displaystyle \frac{1}{10}}\}y^2+{\displaystyle \frac{1}{3}})y^3`$
$`f_{\frac{3}{2}}`$ $`=`$ $`\{[(\{[({\displaystyle \frac{58773}{1114112}}y^2{\displaystyle \frac{77}{5120}})y^2+{\displaystyle \frac{63}{3328}}]y^2{\displaystyle \frac{35}{1408}}\}y^2+{\displaystyle \frac{5}{144}})y^2`$ (38)
$`{\displaystyle \frac{5}{36}}]y^2+{\displaystyle \frac{1}{10}}\}y^5`$
$`f_{\frac{5}{2}}`$ $`=`$ $`[(\{[({\displaystyle \frac{275913}{4456448}}y^2+{\displaystyle \frac{7}{512}})y^2{\displaystyle \frac{7}{416}}]y^2+{\displaystyle \frac{15}{704}}\}y^2{\displaystyle \frac{1}{36}})y^2`$ (39)
$`+{\displaystyle \frac{1}{28}}]y^7`$
### C.2 Small degeneracy ($`\eta <0`$) Cox ; Pichon

Define $`\stackrel{~}{K}_\nu (x)=\mathrm{}^xK_\nu (x)`$ where $`K_\nu `$ is the modified Bessel function, then:

$$F_{\frac{1}{2}}(\eta ,\beta )=\frac{1}{\sqrt{2\beta }}\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^{n1}}{n}\mathrm{}^{n\eta }\stackrel{~}{K}_1\left(\frac{n}{\beta }\right)$$
(40)
$$F_{\frac{3}{2}}(\eta ,\beta )=\frac{1}{\sqrt{2\beta ^3}}\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^{n1}}{n}\mathrm{}^{n\eta }\left[\stackrel{~}{K}_2\left(\frac{n}{\beta }\right)\stackrel{~}{K}_1\left(\frac{n}{\beta }\right)\right]$$
(41)
$$F_{\frac{5}{2}}(\eta ,\beta )=\frac{1}{\sqrt{2\beta ^5}}\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^{n1}}{n}\mathrm{}^{n\eta }\left\{2\left[\stackrel{~}{K}_1\left(\frac{n}{\beta }\right)\stackrel{~}{K}_2\left(\frac{n}{\beta }\right)\right]+\frac{3\beta }{n}\stackrel{~}{K}_2\left(\frac{n}{\beta }\right)\right\}$$
(42)
While for $`\eta >30`$ and $`\mathrm{log}\beta <1`$, the following series expansions by Miralles & Van Riper Mira  are better:

$$F_{\frac{1}{2}}(\eta ,\beta )=\frac{\sqrt{\pi }}{2}\mathrm{}^\eta \left(1+\frac{3}{8}\beta \frac{15}{128}\beta ^2+\frac{105}{1024}\beta ^3\frac{105}{1024}\beta ^4\right)$$
(43)
$$F_{\frac{3}{2}}(\eta ,\beta )=\frac{3\sqrt{\pi }}{4}\mathrm{}^\eta \left(1+\frac{5}{8}\beta \frac{35}{128}\beta ^2\frac{2345}{16384}\beta ^3\right)$$
(44)
$$F_{\frac{5}{2}}(\eta ,\beta )=\frac{\sqrt{\pi }}{2}\mathrm{}^\eta \left(1+\frac{7}{8}\beta \frac{539}{4090}\beta ^2\right)$$
(45)
### C.3 Global fit from EFF Egg ; Pichon

If one changes the parameters of the generalized FD functions from ($`\eta ,\beta `$) to ($`f,g`$) where:

$$\eta =2u+\mathrm{ln}\left(\frac{u1}{u+1}\right)$$
(46)
$$\beta =g/u$$
(47)
and

$$u=\sqrt{1+f}$$
(48)
then the generalized FD functions can be approximated as

$$F_{\frac{1}{2}}(\eta ,\beta )=(\sqrt{2}\beta ^{3/2})^1(\widehat{\rho }+\widehat{u}3\widehat{p})$$
(49)
$$F_{\frac{3}{2}}(\eta ,\beta )=(\sqrt{2}\beta ^{5/2})^1(3\widehat{p}\widehat{u})$$
(50)
$$F_{\frac{5}{2}}(\eta ,\beta )=(\sqrt{2}\beta ^{7/2})^1(2\widehat{u}3\widehat{p})$$
(51)
where

$$\widehat{\rho }=\frac{f}{1+f}g^{3/2}(1+g)^{3/2}\frac{\underset{m=0}{\overset{M}{}}\underset{n=0}{\overset{N}{}}\widehat{\rho }_{mn}f^mg^n}{(1+f)^M(1+g)^N}$$
(52)
$$\widehat{p}=\frac{f}{1+f}g^{5/2}(1+g)^{3/2}\frac{\underset{m=0}{\overset{M}{}}\underset{n=0}{\overset{N}{}}\widehat{P}_{mn}f^mg^n}{(1+f)^M(1+g)^N}$$
(53)
$$\widehat{u}=\frac{f}{1+f}g^{5/2}(1+g)^{3/2}\frac{\underset{m=0}{\overset{M}{}}\underset{n=0}{\overset{N}{}}\widehat{U}_{mn}f^mg^n}{(1+f)^M(1+g)^N}$$
(54)
The values of the coefficients $`\widehat{\rho }_{mn}`$, $`\widehat{P}_{mn}`$ and $`\widehat{U}_{mn}`$ can be found in Egg  with $`M=N=1,2,3`$ or $`4`$.

## Appendix D APPENDIX: Proof of Eq. (9)

We use the series expansion approximation to prove the following expression, for which $`k>1`$.

$`{\displaystyle \frac{\{[F_k(\eta ,\beta )/\eta ]+\beta [F_{k+1}(\eta ,\beta )/\eta ]\}^2}{[F_k(\eta ,\beta )+\beta F_{k+1}(\eta ,\beta )]\{[^2F_k(\eta ,\beta )/\eta ^2]+\beta [^2F_{k+1}(\eta ,\beta )/\eta ^2]\}}}`$
$`\{\begin{array}{cc}1\hfill & \eta \mathrm{}\hfill \\ \frac{k+1}{k}\hfill & \eta \mathrm{},\beta =0,k0\hfill \\ \frac{k+5/2}{k+3/2}\hfill & \eta \mathrm{},\beta 0\hfill \end{array}`$ (58)
### D.1 Small degeneracy, non-relativistic

When $`\eta \mathrm{}`$, $`\beta =0`$, the approximation of the generalized FD functions is given by Cox & Giuli Cox  Eq. (24.111) as:

$$F_k(\mathrm{},0)=F_k(\mathrm{})=\mathrm{}^\eta \mathrm{\Gamma }(k+1)$$
(59)
So we get:

$$F_k(\mathrm{})=\frac{F_k(\mathrm{})}{\eta }=\frac{^2F_k(\mathrm{})}{\eta ^2}$$
(60)
Hence Eq. (D) can be written as:

$$\frac{[F_k(\mathrm{})/\eta ]^2}{F_k(\mathrm{})[^2F_k(\mathrm{})/\eta ^2]}\mathrm{\hspace{0.25em}\hspace{0.25em}1}$$
(61)
### D.2 Small degeneracy, arbitrary relativistic

When $`\eta \mathrm{}`$, $`\beta 0`$, we can find the leading term of the generalized FD functions for the following two extreme cases from Cox  Eq. (24.253) - (24.255) and Eq. (24.267) - (24.269), respectively, as:

$$\{\begin{array}{ccc}F_k(\mathrm{},\beta )=\mathrm{}^\eta g_k(\beta )\hfill & \beta 1\hfill & \\ F_k(\mathrm{},\beta )=\mathrm{}^\eta h_k(\beta )\hfill & \beta 1\hfill & \end{array}$$
(62)
For which we can see that:

$$F_k(\mathrm{},\beta )=\frac{F_k(\mathrm{},\beta )}{\eta }=\frac{^2F_k(\mathrm{},\beta )}{\eta ^2}$$
(63)
Hence Eq. (D) is satisfied.

### D.3 Great degeneracy, non-relativistic

When $`\eta \mathrm{}`$, $`\beta =0`$ we get from Cox  Eq. (24.116) that

$$F_k(\mathrm{},0)=F_k(\mathrm{})=\frac{1}{k+1}\eta ^{k1}$$
(64)
So,

$$\frac{[F_k(\mathrm{})/\eta ]^2}{F_k(\mathrm{})[^2F_k(\mathrm{})/\eta ^2]}\frac{k+1}{k}$$
(65)
### D.4 Great degeneracy, arbitrary relativistic

When $`\eta \mathrm{}`$, $`\beta 0`$, the leading term of the generalized FD functions from Cox  Eq. (24.179) - (24.181) is:

$$F_k(\mathrm{},\beta )=\frac{1}{\sqrt{2}}\eta ^{k+1}\frac{f_k(y)}{(\sqrt{1+y^2}1)^{k+1}}$$
(66)
where $`1+y^2=(1+\eta \beta )^2`$. When $`\eta \mathrm{}`$, $`\beta 0`$, $`y\eta \beta 1`$. So, $`\sqrt{1+y^2}1y`$. From Cox  Eq. (24.202), we can also find

$$f_k(y)\frac{1}{k+3/2}y^{k+3/2}$$
(67)
So,

$$F_k(\mathrm{},\beta )\frac{1}{\sqrt{2}}\eta ^{k+1}\frac{1}{k+3/2}\frac{y^{k+3/2}}{y^{k+1}}=\sqrt{\frac{\beta }{2}}\frac{1}{k+3/2}\eta ^{k+3/2}$$
(68)
Hence we can see that

$$\frac{[F_k(\mathrm{},\beta )/\eta ]^2}{F_k(\mathrm{},\beta )[^2F_k(\mathrm{},\beta )/\eta ^2]}\frac{k+3/2}{k+1/2}$$
(69)
If we look at Eq. (D) we can see that when $`\eta \mathrm{}`$ the $`F_{k+1}(\mathrm{},\beta )`$ terms dominate over the $`F_k(\mathrm{},\beta )`$ terms. As a result, the ratio in Eq. (D) will be $`\frac{k+5/2}{k+3/2}`$.

Table 1: Index of the function or derivative to be evaluated by the program.

| IB | Func. | IB | Func. |
| --- | --- | --- | --- |
| 0 | $`F_k(\eta ,\beta )`$ | 1 | $`\frac{F_k(\eta ,\beta )}{\eta }`$ |
| 2 | $`\frac{F_k(\eta ,\beta )}{\beta }`$ | 3 | $`\frac{^2F_k(\eta ,\beta )}{\eta ^2}`$ |
| 4 | $`\frac{^2F_k(\eta ,\beta )}{\eta \beta }`$ | 5 | $`\frac{^2F_k(\eta ,\beta )}{\beta ^2}`$ |
| 6 | $`\frac{^3F_k(\eta ,\beta )}{\eta ^3}`$ | 7 | $`\frac{^3F_k(\eta ,\beta )}{\eta ^2\beta }`$ |
| 8 | $`\frac{^3F_k(\eta ,\beta )}{\eta \beta ^2}`$ | 9 | $`\frac{^3F_k(\eta ,\beta )}{\beta ^3}`$ |


text:
# Universal Magnetic Fluctuations with a Field Induced Length Scale

## I introduction

There has recently been considerable interest in the fluctuations of a spatially averaged quantity in systems with correlations over a macroscopic length scale $`\xi `$ . The most accessible example from both an experimental and a theoretical point of view is a critical system, where the divergence of the correlation length $`\xi `$ is interrupted by the system size $`L`$. Here, the breakdown of Landau theory and the occurrence of a non-analytic fixed point for the free energy in a renormalization group flow, tells one that the fluctuations will not be Gaussian. Even so, extending the scaling hypothesis to include the macroscopic scale $`L`$, one can deduce that the probability density function (PDF) for order parameter fluctuations should have universal properties. Further, the PDF must be a homogeneous function of $`m`$, $`L`$ and $`\xi `$ of the following form

$`P(m,L,\xi )L^{\beta /\nu }P_L(mL^{\beta /\nu },\xi /L).`$ (1)
Here, adopting the language of a magnetic phase transition, $`\nu `$ and $`\beta `$ are the usual critical exponents relating to the divergence of the correlation length and the singularity in the magnetization, $`m`$. The scaling hypothesis therefore predicts fluctuations of a universal form, independently of system size, for constant ratio $`\xi /L`$.

We study the low temperature phase of the 2D-XY model, defined by the Hamiltonian

$`H=J{\displaystyle \underset{i,j}{}}\mathrm{cos}(\theta _i\theta _j)h{\displaystyle \underset{i}{}}\mathrm{cos}\theta _i.`$ (2)
The exchange interaction and magnetic field are of strength $`J`$ and $`h`$ and the angle $`\theta _i`$ gives the orientation of a classical spin vector of unit length, confined to a plane. We define the magnetization for a single configuration

$`m={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1,N}{}}\mathrm{cos}(\theta _i\overline{\theta }),`$ (3)
where $`\overline{\theta }=\mathrm{tan}^1\left(_i\mathrm{sin}\theta _i/_i\mathrm{cos}\theta _i\right)`$ is the instantaneous magnetization direction.

This is perhaps the simplest nontrivial system in which one can study critical phenomena. At low-temperature and in zero field there is a line of critical points, separated from the high temperature paramagnetic phase by the Kosterlitz-Thouless-Berezinskii phase transition. The physics of this low temperature phase is perfectly captured by a harmonic, or spin wave Hamiltonian. That is, one can, without loss of generality, develop the cosine interaction to order $`(\theta _i\theta _j)^2`$ and neglect the periodicity of $`\theta _i`$. This Hamiltonian is diagonal in reciprocal space and can be solved straightforwardly. As a result, all critical phenomena can be calculated microscopically from Gaussian integration, without the need for either the scaling hypothesis, or the renormalization group. Along the line of critical points the exponents $`\beta `$ and $`\nu `$ are not individually defined, but there ratio is and the system has a single independent exponent $`\eta =2\beta /\nu =T/2\pi J`$.

We have previously been interested in the zero field, or strongly correlated regime where the divergence of $`\xi `$ is completely removed by the system size $`L`$ and $`P_L`$ becomes a function of a single variable $`mL^{\beta /\nu }`$ . We have found that $`P_L`$, when plotted as a function of $`\mu =(mm)/\sigma `$, is a universal function, not only of system size, but also of temperature and therefore of critical exponent $`\eta `$. Here, $`m`$ is the mean and $`\sigma `$ the standard deviation of the distribution. This rather surprising result gives weight to our conjecture  that the critical fluctuations of systems in certain universality classes are captured, at least qualitatively by the fluctuations of the 2D-XY model.

In this paper we generalize our previous results by introducing a second length into the problem with the aid of a magnetic field. The field breaks the symmetry moving the system into an ordered magnetic state with finite correlation length $`\xi `$. However, taking a van Hove type thermodynamic limit, with the ratio $`\xi /L`$ constant, should lead to a family of limit functions, all with divergent correlation length, varying in form from the anisotropic limit (see figure 1) to a Gaussian function, as the ratio $`\xi /L`$ falls to zero. In the next section we develop a starting Hamiltonian that satisfies the requirements of the scaling hypothesis, using a self-consistent Hartree approximation. In section 3 we give theoretical results for the PDF for finite field and compare our results with Monte Carlo simulation. In section 4 we fit the curve with a generalized form of Gumbel’s first asymptote from extremal statistics.

## II Hartree approximation and the hyperscaling relation

Expanding the field energy in small angles, in the same way as the exchange term and Fourier transforming gives a convenient starting Hamiltonian

$`H=H_0+{\displaystyle \frac{J}{2}}{\displaystyle \underset{\stackrel{}{q}0}{}}\left(\gamma _q+{\displaystyle \frac{h}{J}}\right)\varphi _q^2,`$ (4)
where $`\varphi _q=Re[1/\sqrt{N}_i\theta _i\mathrm{exp}(i\stackrel{}{q}.\stackrel{}{r}_i)]`$ is the real part of the Fourier transformed spin variable. $`\gamma _q=42\mathrm{cos}(q_x)2\mathrm{cos}(q_y)`$ and the sum $`_{\stackrel{}{q}0}`$ is over the Brillouin zone for a square lattice with periodic boundaries, with $`\stackrel{}{q}`$ taking on discrete values $`q_x=(2\pi /L)n_x,q_y=(2\pi /L)n_y`$, $`n_x,n_y=0,2\mathrm{}\sqrt{(}N)`$. Here and throughout the paper we have set the nearest neighbour distance on the lattice equal to unity. Expanding $`\gamma _q`$ for small $`q`$ we can write the Green’s function propagator

$`G(q){\displaystyle \frac{1}{q^2+\xi ^2}},\xi =\sqrt{J/h};`$ (5)
the magnetic field indeed introduces a length scale $`\xi `$.

However, this naive starting point needs some development before proceeding with the calculation as, as it stands it does not satisfy the well known hyperscaling relation. To see this, consider the following scaling argument : at the critical temperature but in finite field, the thermally averaged magnetization can be expressed in terms of both $`\xi `$ and $`h`$

$`m\xi ^{\beta /\nu }\left({\displaystyle \frac{h}{J}}\right)^{1/\delta },`$ (6)
where $`\delta `$ is the usual critical exponent for the response in finite field. Putting the expression (5) for $`\xi `$ in (6) leads to a relation between the exponents: $`\delta =2\nu /\beta `$, in disagreement with the hyperscaling relation $`\delta +1=d\nu /\beta `$, which should be valid for the 2D-XY model . The error comes from the development of the field term in small angles. Even at low temperature, when the nearest neighbour differences $`\theta _i\theta _j`$ are small, the deviations of $`\theta _i`$ from the fixed field direction are divergent in the thermodynamic limit. The development of the field term in small angles is therefore invalid. This problem can be dealt with in the low temperature phase, in the absence of vortices, using the Hartree approximation introduced by Pokrovsky and Uimin . Expanding $`\mathrm{cos}(\theta _i)`$ in powers of $`\theta _i^2`$ we make a mean field decoupling

$$\theta _{i}^{}{}_{}{}^{2p}𝒞_p\theta _{i}^{}{}_{}{}^{2p2}\theta _{i}^{}{}_{}{}^{2}$$
(7)
where $`𝒞_p=\frac{(2p)!}{2(2p2)!}`$ is a binomial counting factor. As the underlying Hamiltonian (4) is quadratic we can reduce $`\theta _{i}^{}{}_{}{}^{2p2}`$ using Wick’s theorem and after some resummation we eventually find

$`\mathrm{cos}\theta _i1m{\displaystyle \frac{\theta _i^2}{2}}.`$ (8)
The field term in the Hamiltonian (4) then becomes

$`h{\displaystyle \underset{i}{}}\mathrm{cos}(\theta _i)`$ $`=`$ $`Nh{\displaystyle \frac{1}{2}}h_{eff}(T){\displaystyle \underset{\stackrel{}{q}0}{}}\varphi _q^2`$ (9)
$`h_{eff}(T)`$ $`=`$ $`mh.`$ (10)
Using the scaling relation $`m(h/J)^{\frac{1}{\delta }}`$ the effective field $`h_{eff}h^{\frac{\delta +1}{\delta }}`$ and the above scaling argument correctly yields the hyperscaling scaling relationship defined above. Note however, that the scaling argument is valid in the thermodynamic limit where $`\xi /L1`$ and the influence of the finite system size is negligible. In the crossover region that interests us, with $`\xi /LO(1)`$ one cannot make this substitution and in general one must explicitly work with expression (9).

The point of principle that poses the problem for the hyperscaling relation is that making the substitution $`\mathrm{cos}(\theta _i)1(1/2)\theta _i^2`$ results in an order parameter conjugate to the field, $`m=1(1/2)\theta _i^2=1T/8\pi J\mathrm{log}(N)`$ which diverges with system size, for any finite temperature. In order for the hyperscaling relation to hold, $`m`$ must be a correctly defined intensive variable. For this to be so, the higher order terms in the expansion of $`\mathrm{cos}\theta _i`$ must be retained, at least within the level of approximation shown here. The calculation of Berezinskii  is consistent with this thermodynamic argument.

Still in the absence of vortices, an effective coupling constant $`J_{eff}`$ can be calculated in a similar manner. In zero field we did not need to take account of this, as $`P_L`$ for a quadratic Hamiltonian is independent of temperature throughout the low temperature regime. However, as the correlation length depends on both the field and the coupling constant, we now have to calculate it if we are to have good agreement with numerical data. Expanding $`\mathrm{cos}(\theta _i\theta _j)`$ in powers of the discrete difference operator $`\stackrel{}{}\theta _i`$  and again using the decoupling (7) for $`(\stackrel{}{}\theta _i)^{2p}`$, we arrive at a self-consistent expression for $`J_{eff}(T)`$

$$J_{eff}=J\mathrm{exp}\left(\frac{T}{4J_{eff}}\right).$$
(11)
The Green’s function we finally use for the calculation of the PDF is therefore of the form (5) but with the correlation length given in terms of the self-consistent effective field and coupling constant

$`\xi =\sqrt{{\displaystyle \frac{J_{eff}}{h_{eff}}}}.`$ (12)
The variation of $`\xi `$ with temperature, for fixed field is quite small throughout the range of fields that interests us. In Figure 1 we show $`(\xi (0)/\xi (T))^2`$ as a function of temperature for three different field strengths. Even for $`h/J=0.5`$, there is only a $`10\%`$ variation, up to a temperature $`T/J=0.7`$, above which the Hartree approximation breaks down. As shown below, the temperature dependence of the resulting distribution function is even weaker than that for $`\xi /L`$ and for practical purposes it can be considered as temperature independent.

The parameters $`h_{eff}`$ and $`\xi /L`$ can be found in table 1 for a system of size $`L=32`$ at temperature $`T/J=0.7`$ and for field strength between $`h=0.001`$ and $`h=0.5`$. For finite field the ratio $`\xi /L`$ varies from $`\xi /L=0.99`$ to $`\xi /L=0.043`$.

| $`h`$ | $`h_{eff}`$ | $`\xi /L`$ | $`a`$ | $`\gamma `$ |
| --- | --- | --- | --- | --- |
| $`0.0`$ | $`0.0`$ | $`\mathrm{}`$ | $`1.5807`$ | $`0.89`$ |
| $`0.001`$ | $`0.001`$ | $`0.987`$ | $`1.611`$ | $`0.88`$ |
| $`0.005`$ | $`0.005`$ | $`0.4407`$ | $`1.7416`$ | $`0.844`$ |
| $`0.01`$ | $`0.0101`$ | $`0.3111`$ | $`1.903`$ | $`0.801`$ |
| $`0.05`$ | $`0.0512`$ | $`0.1381`$ | $`3.359`$ | $`0.583`$ |
| $`0.1`$ | $`0.1034`$ | $`0.098`$ | $`5.333`$ | $`0.463`$ |
| $`0.5`$ | $`0.523`$ | $`0.0429`$ | $`21.82`$ | $`0.23`$ |

TABLE 1 Variation of $`h_{eff}`$, $`\xi /L`$, $`a`$ and $`\gamma `$ with field $`h`$ for $`L=32`$

## III The probability density function in finite field

We have previously developed  the following expression for the PDF

$`P_L(\mu )`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}e^{ix\mu }\phi (x)`$ (13)
$`\mathrm{ln}\phi (x)`$ $`=`$ $`ix\sqrt{{\displaystyle \frac{1}{2g_2}}}{\displaystyle \underset{𝐪\mathrm{𝟎}}{}}{\displaystyle \frac{G(𝐪)}{N}}`$ (14)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐪\mathrm{𝟎}}{}}\mathrm{ln}\left[1i\sqrt{{\displaystyle \frac{2}{g_2}}}{\displaystyle \frac{G(𝐪)}{N}}x\right]`$ (15)
where $`g_k=1/N^k_\stackrel{}{q}G(q)^k`$ and $`\mu =(mm)/\sigma `$. The PDF with a finite correlation length is calculated from the same expression by inserting the the modified Green’s function (5). The equivalence of equations (1) and (13) is a result of the hyperscaling relation  and the functional dependence predicted by the scaling hypothesis (1) comes directly from dimensional analysis of equation (13). We note that the calculation can easily be extended to explore the non-Gaussian but non-critical behaviour in all dimension less than four . Summing over the Brillouin zone for a large but finite system and performing a numerical Fourier transform we generate the data shown in figures 2 and 3. Data is shown in figure 2, for $`h/J=0.05`$ ($`\xi /L=0.138`$) and in Figure 3 for $`h/J=0.01,`$ $`0.05`$ and $`0.5`$ corresponding to $`\xi /L=0.311,`$ $`0.138`$ and $`0.043`$. It is compared with results from Monte Carlo simulation for a system of size $`L=32`$. In each case, theoretical and numerical data are shown for three temperatures $`T/J=0.1,\mathrm{\hspace{0.33em}0.4},\mathrm{\hspace{0.33em}0.7}`$. As in the zero field case, which can be considered as the extreme non-Gaussian limit for such a system, the PDFs are characterized by an exponential tail for fluctuations below the mean and a double exponential above the mean . Applying the field reduces the asymmetry and in large field the data approach a Gaussian distribution.

Agreement between the theoretical calculation and the Monte Carlo simulation is generally extremely good, indicating that the Hartree approximations are accurate. In Figure 2 all sets of data collapse, within numerical error, onto a single curve independently of temperature. This is the case for all field values chosen. When plotted on a logarithmic scale, temperature dependence is still not observable, but a difference between theoretical and numerical values can be observed along the exponential tail, for probability densities smaller than $`10^4`$. The discrepancy appears largest for fields around $`h/J0.05`$. This must indicate the limit of the ability of the Hartree approximation in dealing with the fluctuations. For very small field, its effect is small and so errors are negligible, while for larger fields critical fluctuations are smaller and one can imagine that the Hartree approximation becomes quantitatively very accurate. Only in the intermediate field range of $`h/J0.05`$, the combination of these two effects is sufficient to give the small deviation from the Monte Carlo data.

In all cases, both numerical and theoretical, the independence of the results on temperature is quite remarkable and it leads one to suggest that, as in the case of zero field, the distribution is truly temperature independent throughout the range of temperature and system size for which the excitation of vortex pairs can be neglected. This point requires further study, but it is already clear that, from a pragmatic point of view, temperature dependence is not an observable phenomena.

Note that with the definition (3) we explicitly study the longitudinal magnetization, irrespective of its direction in space. This is the quantity that becomes critical at a phase transition in a system with continuous symmetry. The introduction of the magnetic field breaks the orientational symmetry in the thermodynamic limit and the variable conjugate to it is the projection of $`\stackrel{}{m}`$ along the field direction. The fluctuations in these two quantities are different for small field and in the crossover region, as $`h0`$, the latter quantity becomes ill-defined . The two quantities become indistinguishable for $`\xi /L0.1`$.

The skewness $`\gamma =\mu ^3`$, which parameterizes the asymmetry of the curve, varies from $`0.89`$ in the extreme non-Gaussian limit to $`0.23`$ for a ratio $`\xi /L=0.043`$ and indeed goes smoothly to zero for $`\xi /L0`$. Numerical values can also be found in table 1.

The origin of the skewness becomes clear if one considers the contribution made by the normal modes. The magnetization can be written

$`m=11/2N{\displaystyle \underset{q0}{}}\varphi _q^2+\mathrm{}.`$ (16)
To leading order in $`\varphi _q`$ $`m`$ therefore consists of a sum over variables $`m_𝐪=(\varphi _q^2/2N)`$ which, within the spin wave approximation, are statistically independent, with generating function

$`P(m_𝐪)=\sqrt{{\displaystyle \frac{\beta Jq^2N}{4\pi }}}m_𝐪^{1/2}\mathrm{e}^{\beta J_{eff}N(q^2+\xi ^2)m_𝐪}.`$ (17)
In zero field the mean amplitudes $`m_𝐪`$ vary from a microscopic value, $`O(1/N)`$ for modes on the zone edge through to a value of $`O(1)`$ for the long wavelength modes at the zone center and the dispersion in the contributions is divergent in the thermodynamic limit. In two-dimensions the density of states is linear in $`q`$ which is just what is required to engage the entire zone .

Violation of the central limit theorem therefore arises because the individual elements, although statistically independent are not individually negligible. The modes of divergent amplitude near the Brillouin zone center are responsible for the anisotropy although all parts of the zone are required for a detailed reconstruction of $`P(m)`$. Introduction of the length scale $`\xi `$ removes the divergence for $`q0`$ and re-establishes the criterion that the statistically independent elements are individually negligible. In the limit that $`\xi /L0`$ the distribution becomes Gaussian. If the thermodynamic limit is taken while keeping the ratio $`\xi /L`$ constant the amplitudes remains divergent, but the contribution from the zone center becomes progressively less, as the ratio $`\xi /L`$ is reduced and the skewness falls to zero.

## IV Fitting with a generalized Gumbel function

In refs  we have compared the functional form of $`P_L(\mu )`$ with a series of standard expressions. Although none are exact solutions, they all give good fits to the data and therefore offer very useful analytical functions as well as giving some insight into the physical processes responsible for the asymmetric PDF. Here, we only pursue one of these, the generalized Gumbel function for the statistics of extremes

$`\sigma _zP_G(\mu _z)`$ $`=`$ $`w\mathrm{e}^{ab(\mu _zs)a\mathrm{e}^{b(\mu _zs)}},`$ (18)
which gives the $`a^{th}`$ largest or smallest values of a set of $`N`$ random numbers, $`z_i`$, in the limit that $`N\mathrm{}`$. For the smallest values, $`\mu _z=(zz)/\sigma _z`$. The constants $`w`$, $`b`$ and $`s`$ depend on $`a`$ through the three conditions of normalization, $`\mu _z=0`$ and $`\mu _z^2=1`$ and one finds

$`w`$ $`=`$ $`{\displaystyle \frac{a^a\alpha _a}{\mathrm{\Gamma }(a)}}\sigma _z`$ (19)
$`b`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{^2\mathrm{\Gamma }(a)}{a^2}}\left[{\displaystyle \frac{1}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{\mathrm{\Gamma }(a)}{a}}\right]^2}`$ (20)
$`s`$ $`=`$ $`{\displaystyle \frac{1}{b}}\left[\mathrm{log}(a){\displaystyle \frac{1}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{\mathrm{\Gamma }(a)}{a}}\right].`$ (21)
The function therefore has only one parameter, which is calculated by comparing the Fourier transform of (18) with $`\mathrm{\Phi }(x)`$ of equation (13). Using the notation of equation (13) we have, for the Gumbel function

$`\mathrm{ln}\phi _G(x)`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{w\mathrm{\Gamma }(a)}{sa^a}}ix\left(s+{\displaystyle \frac{\mathrm{\Psi }(a)}{b}}{\displaystyle \frac{\mathrm{ln}(a)}{b}}\right)`$ (22)
$``$ $`{\displaystyle \frac{x^2}{2b^2}}\mathrm{\Psi }^{}(a)+i{\displaystyle \frac{x^3}{6b^3}}\mathrm{\Psi }^{\prime \prime }(a)`$ (23)
$`+`$ $`{\displaystyle \frac{x^4}{24b^4}}\mathrm{\Psi }^{\prime \prime \prime }(a)i{\displaystyle \frac{x^5}{120b^5}}\mathrm{\Psi }^{(4)}(a)+\mathrm{}`$ (24)
where $`\mathrm{\Psi }(z)`$ is the digamma function $`\mathrm{\Gamma }^{}(z)/\mathrm{\Gamma }(z)`$. For $`\mathrm{\Phi }`$ we have:

$`\mathrm{ln}\phi (x)={\displaystyle \frac{1}{2}}x^2i{\displaystyle \frac{\sqrt{2}g_3}{3g_2^{3/2}}}x^3`$ (25)
$`+{\displaystyle \frac{g_4}{2g_2^2}}x^4+i{\displaystyle \frac{2\sqrt{2}g_5}{5g_2^{5/2}}}x^5+\mathrm{},`$ (26)
from which it follows that the constant $`a`$ is implicitly given by

$`{\displaystyle \frac{\mathrm{\Psi }^{\prime \prime }(a)}{\mathrm{\Psi }^{}(a)^{3/2}}}=2^{3/2}{\displaystyle \frac{g_3^{3/2}}{g_2}}.`$ (27)
For zero field the solution is $`a\pi /2`$, rather than an integer value, showing $`P(m)`$ is not simply an extreme value distribution. That this solution is an approximation, though be it a good one, can be seen by comparing the ratio of higher order terms in the two expansions (22) and (25). These diverge slowly from unity . Solving equation (27) for finite field gives $`a`$ increasing with $`h`$. The subsequent curves are superimposed in Figures 2 and 3 where one can see that the fitting function reproduces the results of the theoretical calculation to a good approximation. For small field a very small difference in the slopes of the exponential tails can be detected. However, this disappears with increasing $`h`$ and the fitting function can be regarded as an excellent working tool for describing the data. In Figure 4 we illustrate the evolution of the distribution from the anisotropic limit to the uncorrelated Gaussian limit as a function of field using equation (18). The values of $`a(h)`$ are shown in table 1. In terms of extremal statistics, the evolution of $`a(h)`$ means that we are describing the PDF of less and less extreme values, which becomes more and more normal.

For strong field one can solve equation(27) analytically. Evaluating $`g_2`$ and $`g_3`$ using a continuum approximation and using Stirling’s formula, $`\mathrm{ln}\mathrm{\Gamma }(a)a\mathrm{ln}aa`$, one finds

$$a\frac{\pi }{2}\left[1+\left(\frac{L}{2\pi \xi }\right)^2\right].$$
(28)
This simple expression gives $`a`$ to a good approximation even outside the range of $`a`$ values for which Sterling’s formula is valid and reproduces our previous result for $`h=0`$. It also allows one to see that $`a`$ has a contribution coming from fluctuations within a correlated domain and a contribution coming from the fact that the system, with finite $`\xi `$, can be divided into a number, $`N_{eff}=(L/\xi )^2`$, of statistically independent domains.

## V conclusion

In conclusion, we have made a microscopic calculation of the generalized scaling function $`P_L(mL^{\beta /\nu },\xi /L)`$ for order parameter fluctuations near the line of critical points of the low temperature phase of the two-dimensional XY model. A Hartree approximation is used to treat the non-linear corrections to a quadratic Hamiltonian. The approximation is necessary to ensure that the hyperscaling relation between critical exponents is satisfied. We show that the hyperscaling relation is a consequence of the non-linearity necessary to ensure the correct system size dependence of the order parameter, conjugate to the applied magnetic field. This is a requirement of thermodynamics, rather than a general requirement for the observation of non-Gaussian fluctuations for global quantities. Indeed, observation of hyperscaling in non-thermodynamic systems , could be taken as an indication that an equivalent phenomenology exists. For fixed magnetic field the correlation length is modified only slightly by thermal fluctuations and this manifests itself in the function $`P_L`$ which is essentially independent of temperature. The fact that the exponent $`\beta /\nu =4\pi T/J`$ is small and and that $`\delta =8\pi T/J1`$ is large may be important for this observation. More work is required to clarify this point.

Finally, we propose that our fitting parameter $`a(h)`$ could be used as an experimental tool to estimate the correlation length scale in other correlated systems. We have previously made an empirical observation  that the fluctuation of global measures in other correlated systems, both in equilibrium and out of equilibrium are very similar to those of the magnetization of the 2D-XY model in zero field. We have proposed that these observations illustrate at least qualitatively, universal features for correlated systems from different universality classes. This idea, as well as alternative interpretations , could be tested for example, in an enclosed turbulent flow using the experimental set up described in , by varying the ratio of the power injection length scale to the enclosure length scale.

It is a pleasure to thank our collaborators J.-Y. Fortin and J.-F. Pinton for their contributions to this work and L. Berthier, S. Fauve, A. Noullez and Z. Rácz for useful discussions. This work was supported by the PSMN at the École Normale Supérieure de Lyon. M. S. is supported by E.U. contract ERBFMBICT983561.


text:
# On the X-ray fast-time variability of Sco~X-2 (GX~349+2)

## 1 Introduction

Sco~X-2 (GX~349+2) is a neutron star low-mass X-ray binary and Z source. Z sources trace out a characteristic ‘Z’ shaped track in an X-ray colour-colour diagram (CD) or hardness-intensity diagram (HID) consisting of a ‘horizontal’ branch (HB), ‘normal’ branch (NB) and ‘flaring’ branch (FB) (Hasinger & van der Klis hv89 (1989); hereafter HK; van der Klis v95 (1995)). The orientation of each branch in a CD or HID depends upon the choice of energy bands. Mass-accretion rate ($`\dot{M}`$) is thought to increase from the HB, down the NB and up the FB. The other traditional Z sources are: Sco~X-1, GX~17+2, GX~340+0, Cyg~X-2, and GX~5-1.

A class of objects that are related to Z sources are the ‘atoll’ sources. Atoll sources trace out a curved track in a CD, with inferred mass-accretion rate increasing along the track. Atoll sources are thought to have lower magnetic field strengths, and lower mass-accretion rates than Z sources (HK).

The X-ray fast timing properties of Z and atoll sources are different, yet related. The features observed in the power spectra of both classes are normally correlated with the position of the source in a CD or HID. Z sources exhibit very low frequency noise (VLFN), low frequency noise (LFN), high frequency noise (HFN), and quasi-periodic oscillations (QPOs). QPOs with frequencies in the range 15–60 Hz (called HBOs) and 200–1200 Hz (kHz QPOs; these typically occur in pairs) are present in the HB and upper part of the NB. Atoll sources exhibit VLFN, HFN, 1–70 Hz QPOs and kHz QPOs. Atoll source HFN can sometimes be peaked, and resemble broad QPOs (HK and references therein). A recent review of the fast-time variability of Z and atoll sources, with a discussion of various models, has been given by van der Klis (v00 (2000)).

In Z sources, narrow QPOs (called NBOs) with centroid frequencies in the range 5–8 Hz have been observed in the NB. NBOs typically become detectable halfway down the NB, and their properties differ between objects. In Sco~X-1 the frequency of the NBO is $``$5.5–8 Hz in the NB, and jumps to $``$12 Hz as the source enters the FB (the QPO is then referred to as an FBO) (Dieters & van der Klis dv00 (2000)). As the source moves up the FB, the frequency and full width at half maximum (FWHM) increase; eventually, at $``$10 $`\%`$ of the way up the FB, the FBO becomes so broad as to be indistinguishable from the underlying noise (but the relative width never exceeds $``$0.5). In Sco~X-1 the fractional rms increases with increasing $`\dot{M}`$ and is greater at higher energies; and there is no time lag between low and high energy bands (Dieters et al. dvk00 (2000)). The N/FBO seen in GX~17+2 behaves in a similar fashion to Sco~X-1 (Penninx et at. plm90 (1990)). In the other Z sources NBOs have also been observed but are generally less prominent (e.g. Kuulkers et al. kvo94 (1994)).

Sco~X-2 is seemingly very similar to Sco~X-1. The Z tracks have nearly identical shapes (except for the absence of an HB in Sco X-2) and the lightcurves of both objects exhibit strong flaring behaviour (Schultz et al. sht89 (1989); Kuulkers & van der Klis 1995a ). The orbital periods are also similar: 18.9 h for Sco~X-1 (Cowley & Crampton cc75 (1975)); and $``$22 h for Sco~X-2 (Wachter w97 (1997)). One might therefore expect the X-ray fast-time variability in Sco~X-2 to resemble that seen in Sco~X-1; but this is not the case.

EXOSAT observations of Sco~X-2 reported by Ponman et al. (pcs88 (1988)) showed a complete FB. No narrow NBO was found; however broad peaked noise was detected with a centroid frequency and FWHM of about 5 Hz and 10 Hz respectively. The peaked noise was strongest in the lower FB, with a fractional rms of 3.2$`\pm `$0.2 $`\%`$. There was no significant change in the frequency or width of the peaked noise, and there was no time lag between the low and high energy bands. We note that the authors divided their data into four segments on the basis of intensity, and not according to position in the Z track. By using that method they could not differentiate between the NB and the lower FB.

A small set of Rossi X-ray Timing Explorer observations were reported by Kuulkers & van der Klis (kv98 (1998); hereafter KK), when Sco~X-2 was near the NB/FB vertex. In the lower NB they found peaked noise with a centroid frequency, FWHM and fractional rms of 9.4$`\pm `$0.5 Hz, 16$`\pm `$2 Hz and 2.7$`\pm `$0.1 $`\%`$ respectively. In the lower region of the FB they measured these values as 5.8$`\pm `$0.2 Hz, 11.0$`\pm `$0.7 Hz and 4.2$`\pm `$0.1 $`\%`$. There was a small but significant decrease in frequency when the source moved through the NB/FB vertex. They commented that the peaked noise resembled the peaked HFN sometimes seen in atoll sources. They found that the peaked noise fractional rms was stronger at higher energies.

Even though there have been clear indications that Sco~X-2 exhibits marked differences to Sco~X-1, it has remained largely ignored. This is perhaps due to the fact that it looks so similar to Sco~X-1 in a CD and HID.

We present results obtained from an analysis of archived data from the Ginga satellite; these data have not previously been presented. This work is the first detailed analysis of the power spectral properties of Sco~X-2 as a function of position in the Z track. This analysis method is extremely useful because it allows us to unambiguously distinguish between the NB and lower FB, and is helpful when making comparisons between different objects. We further highlight the differences between Sco~X-2 and Sco~X-1, and conclude that the observed peaked noise is not the same as atoll source peaked HFN, as had been suggested. We compare our results with the power spectral features seen in other Z sources.

## 2 Observations

The Ginga Large Area Counter (LAC) (Makino & ASTRO-C Team m87 (1987); Turner et al. ttp89 (1989)) observed Sco~X-2 between 5 and 7 March 1989. We retrieved the data from the Leicester Database and Archive Service111see http://ledas-www.star.le.ac.uk/ginga/. In Table 1 we present a log of observations.

The LAC consisted of eight detectors, and was sensitive to X-rays in the range 1–37 keV. Several modes were available, with differing time and energy resolutions (higher time resolution sacrifices energy resolution). The observations discussed here were done in the so-called MPC3 and PC modes. For the MPC3 mode the coarse gain setting was high, providing an energy range 1–19 keV over twelve channels. The time resolution was 7.8 ms and 62.5 ms for high (MPC3-H) and medium (MPC3-M) bit rates respectively. In the PC mode, the LAC was divided into two groups of four detectors each. The coarse gain setting was low, giving an energy range 1.5–24.3 keV for one detector group, and 1.2–17.9 keV for the other. Each detector group had two energy channels, so the PC mode observations had four overlapping energy bands. In the high bit rate mode (PC-H), the time resolutions were 0.98 ms and 1.95 ms for the lower and higher energy bands respectively. Our 0.98 ms data were rebinned with a time resolution of 1.95 ms, to obtain lightcurves for the full energy range.

During the observations in section C (see Table 1) the high voltages in detectors 5 and 6 were accidently changed. The energy range for the MPC3 observations stated above is therefore only approximate for section C.

## 3 Analysis

The counting rates were corrected for detector deadtime and collimator response. During the observations in sections A-C, the collimator response was $``$0.96. In section D the response was $``$0.73. The data were not corrected for collimator reflection. The effect of collimator reflection was to overestimate counting rates below $``$6 keV by a few percent (Turner et al. ttp89 (1989); see also Hasinger et al. hve90 (1990); Wijnands et al. wvk97 (1997)).

Hardness-intensity diagrams were constructed using 64 s averages. For the MPC3 data we defined the hardness as the counting rate ratio between 8.2–16.7 keV and 5.8–8.2 keV. The intensity was defined as the counting rate in the range 1.1–16.7 keV. For the PC data we defined the hardness as the counting rate ratio between 5.7–17.9 keV and 1.2–5.7 keV. The intensity was defined as the counting rate, from one detector group only (i.e. four detectors), in the range 1.2–17.9 keV.

Timing analysis was performed on the raw (i.e. not corrected for deadtime and collimator response) data. Power spectra were calculated from 16 s intervals in several energy bands, and then averaged according to either corrected counting rate or rank number.

Rank number was introduced by Hasinger et al. hve90 (1990) (see also Lewin et al. llt92 (1992); Hertz et al. hvw92 (1992); Dieters & van der Klis dv00 (2000)) and is a one-dimensional measure of the position of the source in a CD or HID. Rank number increases smoothly with increasing inferred mass-accretion rate. To measure the rank numbers, a spline was fitted through the Z track, and each point was projected onto it. Rank number was then determined by measuring the distance of each point along the spline. We defined the NB/FB vertex as rank number 1, and the top of the observed FB as rank number 2.

We calculated the expected deadtime-affected white-noise levels and then subtracted the white-noise from each average power spectrum. The procedure to calculate white-noise levels in sum channel power spectra is well known, and is described by van der Klis (v89 (1989)) and Mitsuda & Dotani (md89 (1989)). The full procedure to calculate white-noise levels in energy resolved power spectra is not explicity stated in those references; we therefore describe our procedure in the Appendix.

After the deadtime-affected white-noise levels had been subtracted, the spectra were normalised to fractional-rms-squared per Hz and functional $`\chi ^2`$ fits were performed. Very low frequency noise (VLFN) was modelled using a power law $`A\nu ^\alpha `$ and integrated over the range 0.01–1 Hz. Peaked noise was modelled using a Lorentzian. All power spectra were initially fitted with only a power law. The significance of adding a Lorentzian component to the fit was then estimated via an F-test. If the inclusion of a peaked noise component was not significant ($`<3\sigma `$) we measured 90 $`\%`$ confidence upper limits by using $`\mathrm{\Delta }\chi ^2=2.71`$. Errors on all parameters were measured using $`\mathrm{\Delta }\chi ^2=1`$. When a fit was not sensitive to the centroid or FWHM of the Lorentzian (or was not significantly improved by its inclusion) we fixed those parameters at values measured from other power spectra (the actual values are shown in Table 2). We fitted each PC mode power spectrum both with and without an extra free constant component. The reasons for doing this are discussed in Sect. 4.3.

All fractional rms amplitudes were corrected for channel cross-talk (Lewin et al. llt92 (1992)). We assumed variations at all energies were in phase. Finally, the fractional rms values were multiplied by a binning correction factor (van der Klis v89 (1989)).

## 4 Results

### 4.1 Hardness-intensity diagrams

Fig. 1 shows a HID for all section B data. The small, approximately vertical NB is clearly visible. Sco~X-2 is known to flare to intensities roughly twice the level of the persistent emission. The highest intensity in Fig. 1 is a factor of 2.2 greater than the lowest intensity: we conclude therefore, that we have observed a fully developed flaring branch. The results from section C are compatible with those from section B; however due to the uncertain energy ranges of those data, we will not discuss them any further.

A HID of sections A and D is shown in Fig. 2. The lowest and highest intensities in Fig. 2 differ by a factor of 2; therefore (in comparison with Fig. 1), nearly all of the FB was traced out. No clear NB/FB vertex is apparent (but we were able to estimate where it might lie; see Sect. 4.3). When we constructed a HID from MPC3 mode data using energy bands similar to those available in PC mode, the presence of the NB was much less pronounced; the NB became more like a ‘blob’, rather than a clear branch. Therefore we caution that in Fig. 2 we cannot confidently distinguish between data points from the lower NB (if there are any) and those actually in the lower FB. However, as will be seen in Sect. 4.3, our main interest (with the PC data) is with the top half of the FB, so this ambiguity is not a major problem.

### 4.2 Timing analysis of section B MPC3 mode data

During the observations in MPC3-H mode, the source was found only in the lower FB. We divided the data into two segments with mean ranks 1.09$`\pm `$0.02 and 1.26$`\pm `$0.05. An example of a power spectrum (1.1–16.7 keV) is shown in Fig. 3, with mean rank number 1.09. The results of fitting the MPC3-H data (1.1–16.7 keV) are shown in Table 2. The FWHM of the best-fit Lorentzian is greater than half its centroid frequency, so the feature cannot be formally described as a QPO. The MPC3-H data show that there is no significant change in the centroid frequency or FWHM of the peaked noise, as the source moves up the FB (the decrease in the centroid frequency between rank numbers 1.09 and 1.26 is significant only at the $``$$`\sigma `$ level). The fractional rms amplitude of the peaked noise is also consistent with remaining constant.

We investigated the energy dependence of the peaked noise fractional rms, by calculating energy resolved power spectra from the MPC3-H data. We used four energy bands: 1.1–3.4–5.8–8.2–16.7 keV. All power spectra from a particular band were averaged together, giving a mean power spectrum from the lower FB; the mean rank number was 1.18$`\pm `$0.05. We determined the centroid frequency and FWHM of the peaked noise feature from an average power spectrum from the 1.1–16.7 keV band: they were 6.7$`\pm `$0.4 Hz and 8.3$`{}_{1.1}{}^{}{}_{}{}^{+1.2}`$ Hz respectively. The centroid and width were fixed at those values when fitting the power spectra from the other four energy bands. The energy dependence of peaked noise fractional rms is shown in Table 3. The peaked noise is clearly stronger at higher energies and is consistent with being constant above $``$6 keV.

The MPC3-M data from the track shown in Fig. 1 were divided into eight segments. Power spectra were averaged on the basis of rank number. An example (1.1–16.7 keV) is shown in Fig. 4, with mean rank number 1.09. Table 2 shows the results of fitting the MPC3-M power spectra (1.1–16.7 keV) from section B. It was only in the lower FB that we could meaningfully constrain the centroid and width of the peaked noise, which were $``$4 Hz and $``$6 Hz respectively. There is no significant change in those values between mean ranks 1.09 and 1.27 (Table 2). The centroid frequency is lower in the MPC3-M data than in the MPC3-H data. This discrepancy is probably due to a lack of information above 8 Hz in the MPC3-M data.

The dependence of peaked noise and VLFN fractional rms on rank number is shown in Figs. 5 and 6 respectively. The peaked noise fractional rms is greatest in the lower FB and becomes undetectable (upper limit 1.3 $`\%`$) by about halfway up the FB. Of note is the fact that it *is* actually detectable up to a mean rank of 1.41. (This is important because N/FBO typically become undetectable by $``$10 $`\%`$ of the way up the FB.) The drop in power as the source moves further up the FB or into the NB may be due to the feature moving above the Nyquist frequency of these data; however we can exclude this effect with respect to movement up the FB (see Sect. 4.3). The VLFN fractional rms is seen to increase as the source moves up the FB, as expected for Z sources (HK). The slope of the power law appears to be uncorrelated with rank number.

### 4.3 Timing analysis of PC mode data

The NB/FB vertex was not unambiguously observed, so we did not calculate rank numbers for each individual interval; averaging was done on the basis of intensity. The track was divided into four segments; the intensity ranges were: 2800–3200, 3200–3900, 4500–5200 and 5200–5900 count/s (four detectors only). The mean intensities of those four segments were: 3041$`\pm `$6, 3458$`\pm `$14, 4962$`\pm `$25 and 5492$`\pm `$31 counts/s. We could, however, estimate rank numbers by comparing Fig. 2 with Fig. 1 (assuming the intensities are derived from the same energy band). As stated in Sect. 4.1, the track shown in Fig. 2 is not quite a complete FB. The lowest intensity observed in Fig. 2 is $``$2900 counts/s; this can be taken as the upper limit of the intensity of the NB/FB vertex. Alternatively, if we assume the highest intensity present in Fig. 2 is actually the top of the FB, then the unseen NB/FB vertex can be placed at $``$2780 counts/s; this can be considered to be the lower limit. By considering these two extremes we calculated the approximate rank numbers of the four segments as: 1.06, 1.20, 1.68, and 1.80 ($`\pm `$0.03). We stress that these rank numbers are only a rough guide, and that the lower FB in Fig. 2 may infact contain data from the NB.

An example power spectrum is shown in Fig. 7 (1.2–17.9 keV), with mean intensity 3041 counts/s (mean rank number $``$1.06). The results from fitting the PC power spectra (1.2–17.9 keV) are shown in Table 4, where we present the fits from models both with and without a constant. The results from fits with a constant are shown also in Figs. 5 and 6. The most important thing to note is the nondetection of peaked noise in the upper FB (upper limit 1.6 $`\%`$). The PC data Nyquist frequency is 256 Hz; therefore the nondetection of peaked noise in these data is not due only to the feature moving to higher frequencies. From this, we conclude that the nondetection of peaked noise in the MPC3-M data (Table 2) is also a real weakening of the peak.

The reduced $`\chi ^2`$s of the fits without a constant in Table 4 are formally acceptable; however visual inspection of the power spectra revealed excess power above about 100 Hz (e.g. see Fig. 7), and the best-fit VLFN indices found in the two lowest intensity segments are unusually flat. Given these indications, we attempted to account for a possible high frequency feature. We added a cut-off power law to our model and attempted to fit the spectrum from the lowest intensity segment. We could not constrain the index or cut-off frequency of this component. When we fitted a Lorentzian (instead of the cut-off power law) to the same spectrum we could constrain neither the width nor centroid frequency. The fits were improved with the addition of the extra component, but were not sensitive to the exact form of that component. We fitted the spectrum with a constant (instead of a cut-off power law or Lorentzian) and the fit also improved. None of the best-fit values of the peaked noise or VLFN were sensitive to the type of model that we used to describe the excess power; for simplicity, we therefore chose to use a constant.

In the three highest intensity segments, the inclusion of a constant significantly improved the fits (determined via an F-test for the inclusion of the extra component). Furthermore, in the segment with mean intensity 3458 counts/s, fitting without a constant resulted only in determining an upper limit on the peaked noise of $`<`$3.1 $`\%`$ (Table 4). This upper limit is rather large, and a visual inspection of the power spectrum showed peaked noise to be present. Therefore, it seems justified to include an extra component in the fits, without needing (nor being able) to describe the actual form of the possible higher frequency feature. We note that: HFN was reported in Sco~X-2, with cut-off frequencies of 66 Hz and 41 Hz in the NB and FB respectively (HK); also a QPO with a frequency of $``$125 Hz and FWHM of $``$30 Hz was detected by Penninx et al. (plm90 (1990)), in the FB of GX~17+2.

The fact that we cannot determine the form of the higher frequency feature is not a major problem. The main result we obtained from the PC data is this: the peaked noise feature experienced a real weakening in the upper FB. This result is not sensitive to the presence or otherwise of a higher frequency feature.

## 5 Discussion

### 5.1 Observed properties of Sco~X-2

We were able to trace the power spectral features of Sco~X-2 as a function of position in the Z track, in the energy range $``$1–17 keV. At no time did we detect a narrow N/FBO as found in Sco~X-1. Broad peaked noise was present at $``$4–7 Hz, with FWHM $``$6–12 Hz. We will subsequently refer to this feature as flaring branch noise (FBN). FBN was detected in the lower NB and reached a maximum at rank number 1.09 ($``$10 $`\%`$ of the way up the FB). At its maximum, the fractional rms was $``$$`\%`$. It then became weaker and was undetectable by rank number 1.5 (the middle of the FB; upper limit 1.3 $`\%`$). We found no significant change (more than a few Hz) in the centroid frequency or FWHM as the source moved up the FB. The FBN was stronger at higher energies.

The behaviour we have observed is very similar to that reported by Ponman et at. (pcs88 (1988)). In addition, we have been able to unambiguously separate NB data from lower FB data. The increase in the fractional rms that accompanied movement from the NB through to the lower FB, was also seen in the RXTE data of KK. We have not been able to confirm the decrease that they observed in the centroid frequency and FWHM, which accompanied that transition.

### 5.2 Comparison with other power spectral features

The fast-time variability of Sco~X-2 is different to that seen in the other Z sources. Most obvious is the low coherence of FBN when compared with N/FBOs. The FBN has a relative width of $``$1–2, compared with the N/FBO in Sco~X-1 which has a relative width of less than $``$0.5 (Dieters & van der Klis dv00 (2000)). Also in Sco~X-1, the frequency and FWHM of the NBO increase with increasing inferred mass accretion rate, and there is a sudden jump in frequency at the NB/FB vertex. In Sco~X-2 on the other hand, there is a decrease in the FBN centroid frequency and FWHM with movement through the vertex (KK). In Sco~X-1 the NBO strength increases with increasing rank number; the NBO becomes indistinguishable from the underlying noise (by $``$10 $`\%`$ of the way up the FB) because the FWHM increases. In comparison, FBN is strongest at about 10 $`\%`$ of the way up the FB, and becomes undetectable because of a subsequent decrease in the fractional rms.

Even given the marked differences, there are some similarities between FBN and the N/FBO seen in Sco~X-1: the centroid frequencies are similar, both the N/FBO and FBN are hard, and there is no time lag between the soft and hard energies in either phenomenon. In Sco~X-1, the NB/FB vertex is a critical position in the Z track. In Sco~X-2 there are hints that this may also be true, although in a different sense (i.e. frequency *jump* for NBOs compared with a frequency *fall* for FBN). FBN and Z source N/FBOs may still be related, as already suggested by KK.

In the scheme of the unified model, radiation pressure drives matter up out of the disc, forming an inner disc corona, resulting in a radial inflow (Lamb l91 (1991); Psaltis et al. plm95 (1995)). When $`\dot{M}`$ rises to within $``$20 $`\%`$ of the Eddington limit, radial oscillations in the optical depth of the inflowing material produce a rocking of the X-ray spectrum; this gives rise to NBOs (Fortner et al. flm89 (1989); Miller & Park mp95 (1995)). This model may explain the evolution of the N/FBO seen in Sco~X-1, but it does not do well in explaining FBN. The FBN centroid frequency does not increase as predicted by the model. Furthermore, we have detected it $``$40 $`\%`$ of the way up the FB; oscillations are predicted to be suppressed at such high accretion rates.

QPOs at super-Eddington accretion rates have been observed before: Kuulkers & van der Klis (1995b ) observed a $``$26 Hz QPO in the upper FB of Cyg~X-2. The QPO was observed during an intensity dip in the FB, and it was stronger at lower energies. They suggested that it is produced by oscillations in an inner disc torus, and it is observable because Cyg~X-2 is viewed at a high inclination (see Sect. 5.3). Sco~X-2, on the other hand, is thought to be viewed at a low inclination (see Sect. 5.3); therefore it is unlikely that FBN is related to the 26 Hz QPO.

KK pointed out that the FBN seen in Sco~X-2 looks similar to peaked HFN sometimes seen in atoll sources. We observed a *decrease* in the FBN fractional rms as the source moved up the NB, which is opposite to what would be expected if the FBN was infact atoll source HFN. The low Nyquist frequency is no problem here as the break frequency of HFN decreases with decreasing inferred mass-accretion rate (Wijnands & van der Klis wv99 (1999)). The phenomenology of FBN is thus significantly different to atoll source HFN; therefore we propose that FBN is unrelated to atoll source HFN.

### 5.3 Two types of Z sources: a clue to the origin of FBN?

In order to understand why Sco~X-2 behaves so differently, it is worthwhile to investigate the general relationships that exist between the Z sources. On the basis of their Z track morphology and fast-time variability, they can be divided into two groups: Cyg~X-2, GX~5-1, GX~340+0 (Cyg-like objects); and Sco~X-1, Sco~X-2, GX~17+2 (Sco-like objects) (e.g. Kuulkers et al. kvo97 (1997)). The Cyg-like sources have much larger HBs and smaller FBs than the Sco-like sources. The HBs in the Sco-like objects are slanted, while the HBs in the Cyg-like sources are horizontal. Most importantly, in the Cyg-like objects the intensity is seen to fall in the middle and upper FB, while the Sco-like objects really do flare in the FB.

At high accretion rates an inner disc torus may be produced, which in high inclination sources can partly obscure the X-ray emitting region. The Cyg-like objects have thus been interpreted as being viewed at higher inclinations, and the Sco-like objects at lower inclinations (Kuulkers et al. kvo97 (1997)). In support of this, optical observations place Cyg~X-2 at an inclination of 62.5$`\pm `$4$`\mathrm{°}`$ (Orosz & Kuulkers ok99 (1999)), Sco~X-1 at 15$`\mathrm{°}`$–40$`\mathrm{°}`$ (Crampton et al. cch76 (1976)), and also suggest that Sco~X-2 is viewed at a low inclination (Wachter & Margon wm96 (1996)). Psaltis et al. (plm95 (1995)) have modelled the X-ray spectra of Z sources. They suggested that the Sco-like objects have weaker magnetic fields than the Cyg-like objects.

The most obvious difference between Sco~X-2 and the other Z sources is the absence of an HB, so it is tempting to think that either: Sco~X-2 has a weaker magnetic field than Sco~X-1, while also (perhaps incidentally) being viewed at a low inclination; or simply that Sco~X-2 never reaches a low enough mass-accretion rate to trace out an HB.

A close investigation of the timing properties of Cyg~X-2 may provide hints as to the origin of FBN. In Cyg~X-2 long-term variations in the mean X-ray intensity have been discovered (e.g. Paul et al. pkm00 (2000)). Ginga observations of Cyg~X-2 showed that the Z track was more Sco-like at medium level intensities, and more Cyg-like at high level intensities (Wijnands et al. wvk97 (1997)). RXTE observations of Cyg~X-2 at low intensities revealed no NBO (upper limit $``$$`\%`$) in the NB; however a peaked noise feature was detected in the lower FB, with a frequency and fractional rms of 6–7 Hz and $``$$`\%`$ respectively (Kuulkers et al. kwv99 (1999)). This is (perhaps only superficially) similar to the FBN seen in Sco~X-2.

It is vital to have more observations of Cyg~X-2 at low intensities so that we can better understand the properties of the peaked noise feature. Only then can we make any conclusions about the relationship between FBN and the Cyg~X-2 peaked noise. But if the two features are the same phenomenon, we must ask: what are the physical characteristics of Cyg~X-2 when it is at low level intensities? The low level may be due to obscuration by the outer accretion disc, possibly caused by a precessing disc (e.g. Kuulkers et al. kvv96 (1996)). Wijnands & van der Klis (wv01 (2001)) have argued against the precessing disc model on the basis of differing kHz QPO properties at different intensity levels. Kuulkers et al. (kwv99 (1999)) showed that the relative VLFN properties at different intensity levels were also inconsistent with obscuration. Instead, they preferred a model in which the overall intensity level is sensitive to the properties of the inner mass flow. Whatever the case, if FBN and Cyg~X-2 peaked noise are produced by the same process, then that process will be largely insensitive to orbital inclination. In this way, Sco~X-2 may prove to be a useful tool in discerning between models trying to explain the behaviour of Cyg~X-2, and ultimately the other Z sources.

## 6 Conclusions

We have carried out the first detailed study of the X-ray fast-time variability of the Z source Sco~X-2, as a function of position in the Z track. We found broad peaked noise with centroid frequency and FWHM in the range $``$4–7 Hz and $``$6–12 Hz respectively. We find it difficult to attribute FBN to the same processes which are able to explain typical Z source N/FBOs. We also conclude that it is not a manifestation of atoll source HFN. FBN resembles the FB peaked noise seen in the FB of Cyg~X-2 at low overall intensities. If these are later confirmed to be the same phenomena, then the process responsible for FBN must be largely insensitive to orbital inclination.

###### Acknowledgements.

This research has made use of data obtained from the Leicester Database and Archive Service at the Department of Physics and Astronomy, Leicester University, UK. The reduction of raw data was carried out with support from the Cooperative Research Centre for Advanced Computational Systems established under the Australian Government’s Cooperative Research Centres Program. PMO wishes to thank the Space Research Organization Netherlands for its hospitality, and is currently being supported by the Australian Postgraduate Award scheme.

## Appendix: Calculation of deadtime-affected white-noise levels

In the absence of deadtime, the white-noise level in a Leahy normalised power spectrum, due to Poisson noise in the data, is equal to 2 (Leahy et al. lde83 (1983)). The effect of deadtime is to decrease the variance of the data and thus lower the white-noise level. The deadtime-affected white-noise level for the sum channel $`P_{\nu \mathrm{s}}`$, can be calculated for each frequency $`\nu `$ in the power spectrum

$$\begin{array}{c}<P_{\nu \mathrm{s}}>=2(1\mu _\mathrm{s}\tau _{\mathrm{dead}})^2\times \hfill \\ \hfill \left(1+2\left(\frac{\mu _\mathrm{s}\tau _{\mathrm{dead}}}{1\mu _\mathrm{s}\tau _{\mathrm{dead}}}\right)\left(\frac{\tau _{\mathrm{dead}}}{\tau _{\mathrm{res}}}\right)\mathrm{sin}^2\left(\frac{\pi \nu }{2\nu _{\mathrm{Nyquist}}}\right)\right)\end{array}$$
(1)
where $`\mu _\mathrm{s}`$ is the observed counting rate per detector, $`\tau _{\mathrm{res}}`$ is the time resolution of the lightcurve, $`\tau _{\mathrm{dead}}`$ is the deadtime per event, and $`\nu _{\mathrm{Nyquist}}`$ is the Nyquist frequency of each transform (Weisskopf w85 (1985)).

The deadtime-affected white-noise level for any particular energy band $`P_{\nu \mathrm{b}}`$, can be calculated as a function of $`P_{\nu \mathrm{s}}`$. For an arbitrary time interval, the ratio of the total number of counts $`k`$ in the particular energy band, to the total number of counts $`n`$ in the sum channel is denoted by $`r`$. The mean variance of the data is denoted as $`\sigma _\mathrm{s}^2`$ and $`\sigma _\mathrm{b}^2`$ for the sum channel and the particular energy band respectively. These are related by (Mitsuda & Dotani md89 (1989))

$$\sigma _\mathrm{b}^2=r^2\sigma _\mathrm{s}^2+k\left(1r\right)$$
(2)
Each transform is calculated from a time series $`x_k`$ ($`k=0,\mathrm{},N1`$) consisting of $`N`$ arbitrary time intervals. By employing Parseval’s theorem, we can relate the total variance $`Var\left(x_k\right)`$ (or the mean variance $`\sigma ^2`$) of the time series, and the Fourier amplitude $`a`$

$$Var\left(x_k\right)=N\sigma ^2=\frac{N1}{N}|a|^2$$
(3)
Here, we have considered only situations where the dependence of Fourier amplitude on frequency is negligible (i.e. white-noise). The factor of N-1 originates from the fact that the sum used in the expression for the total variance does not include a term for a frequency of zero.

The mean variance $`\sigma ^2`$ and Leahy power $`P`$ can then be related via

$$\sigma ^2=P\frac{N_{\mathrm{ph}}}{2}\left(\frac{N1}{N^2}\right)$$
(4)
$`N_{\mathrm{ph}}`$ is equal to $`Nn`$ for the sum channel or $`Nk`$ for the particular energy band. We can express $`P_\mathrm{b}`$ as a function of $`P_\mathrm{s}`$ by substituting Eq. 4 into Eq. 2

$$P_\mathrm{b}\frac{k}{2}\left(\frac{N1}{N}\right)=P_\mathrm{s}\frac{k^2}{n^2}\frac{n}{2}\left(\frac{N1}{N}\right)+k\left(1r\right)$$
(5)
or

$$P_\mathrm{b}=2\left(\frac{N}{N1}\right)\frac{k}{n}\left(2\left(\frac{N}{N1}\right)P_\mathrm{s}\right)$$
(6)
which in the case of large $`N`$ reduces to

$$P_\mathrm{b}=2\frac{k}{n}\left(2P_\mathrm{s}\right)$$
(7)