text:
# References

I. Introduction

On September 8, 2002, a conjunction of quasar J0842+1835 and Jupiter took place. This event was used to measure the Shapiro time delay of the quasar signal due to the gravity of Jupiter.

Many years ago, I. I. Shapiro proposed one of the classic tests of general relativity in which radio signals are bounced off an inner planet during a superior conjunction with the Sun. The effect of the Sun’s gravity is to create a delay in the time required for the radio waves to return to Earth. In subsequent years, measurements performed using Mercury confirmed Einstein’s theory, and the PPN parameter $`\gamma `$ was measured to be its expected value of $`1`$ to within 10%.

Because Jupiter’s gravity is weaker than the Sun’s, the QSO J0842+1835 measurement required remarkable accuracy: $`10^{12}`$ seconds. This was achieved using very long baseline interferometry. Motivation for undertaking this experiment stems from a proposal that it can be used to measure the speed of gravity $`c_g`$. The idea of testing whether $`c_g`$ equals the speed of light $`c`$, as should be the case in general relativity, has attracted considerable attention both in the astrophysics community and in the media. The measurement yielded $`c_g/c=1.06\pm 0.21`$ and was hailed as a confirmation of Einstein’s general theory of relativity.

The purpose of this Letter is to point out an error in the theoretical formula used to analyze the Jupiter/quasar experiment and to provide the correct result. In reference, a $`v/c_g`$ correction to the Shapiro time delay in the Jupiter/quasar experiment is found to be proportional to $`1/\theta ^2`$, where $`\theta `$ is the angle between the quasar and Jupiter. Since $`\theta `$ is small, an enhancement occurs thereby making the measurement feasible. However, using a simple method, this Letter computes the $`v/c`$ corrections and finds no such term. The discrepancy between the formula of the current work and the one used in the experiment is understood: The angle $`\theta `$ in the latter was actually not the observable one but an artificially defined angle.

Our notation conforms to that of references  and , which are henceforth indicated as A and B: Quasar J0842+1835 is located in the direction of the unit vector $`\stackrel{}{K}`$. See Figure 1. Radiation for the quasar arrives at two observational points $`1`$ and $`2`$ on Earth, which are separated from one another by the distance $`\stackrel{}{B}`$. The impact parameters for each of these two points is respectively denoted by $`\stackrel{}{\xi }_1=\xi _1\stackrel{}{n}`$ and $`\stackrel{}{\xi }_2=\xi _2\stackrel{}{n}`$. Here, $`\stackrel{}{n}`$ is a unit vector perpendicular to $`\stackrel{}{K}`$ going from Jupiter to the closest approach of the electromagnetic radiation of the quasar. Since the difference of the impact parameters is small compared to either impact parameter, we use $`\xi `$ to denote the value of either when a distinction is not important. The velocity of Jupiter is indicated as $`\stackrel{}{v}_J`$, and the Earth-Jupiter distance is denoted by $`R_{EJ}`$.

We are interested in the most significant corrections to the Shapiro time delay for the Jupiter/quasar experiment. Therefore, we neglect terms proportional to the product of two of, or the square of, any of the following small, dimensionless quantities: $`\frac{G_NM_J}{\xi c^2}6\times 10^9`$, $`\frac{v_J}{c}4.5\times 10^5`$, $`\frac{B}{\xi }0.006`$ and $`\theta _{obs}=\frac{\xi }{R_{EJ}}0.001`$ (which is the angle that an astronomer observes between Jupiter and the quasar). Here, $`G_N`$ is Newton’s constant and $`M_J`$ is the mass of Jupiter.

II. The $`v_J/c`$ Corrections

If $`\mathrm{\Delta }t_1`$ and $`\mathrm{\Delta }t_2`$ denote the Shapiro time delays at the points $`1`$ and $`2`$, then the quantity of interest is the difference $`\mathrm{\Delta }(t_1,t_2)=\mathrm{\Delta }t_2\mathrm{\Delta }t_1`$:

$$t_2t_1=|\stackrel{}{x}_2(t_2)\stackrel{}{x}_0|/c|\stackrel{}{x}_1(t_1)\stackrel{}{x}_0|/c+\mathrm{\Delta }(t_1,t_2).$$
(1)
Here, $`t_1`$ and $`t_2`$ are respectively the times at which the signals are measured at the two points $`\stackrel{}{x}_1(t_1)`$ and $`\stackrel{}{x}_2(t_2)`$ on Earth, $`\stackrel{}{x}_0`$ is the position of the quasar, and $`|\stackrel{}{x}_2(t_2)\stackrel{}{x}_0|/c|\stackrel{}{x}_1(t_1)\stackrel{}{x}_0|/c`$ is the time difference that occurs when gravitational effects are absent.

If $`\stackrel{}{B}=\stackrel{}{x}_2(t_2)\stackrel{}{x}_1(t_1)`$ and $`\stackrel{}{n}`$ are oppositely oriented, or more precisely $`\stackrel{}{B}\stackrel{}{n}<0`$, then $`\xi _1>\xi _2`$ and $`\mathrm{\Delta }(t_1,t_2)`$ is positive because the electromagnetic radiation that arrives at $`2`$ undergoes more time delay because it passes closer to Jupiter. This is the case illustrated in Figure 1.

If Jupiter were not moving, which is the static situation, then the Shapiro time delay for a single wave is

$$\mathrm{\Delta }t=\frac{2G_NM_J}{c^3}\left(1+\mathrm{ln}\left(\frac{4R_{JQ}R_{EJ}}{\xi ^2}\right)\right),$$
(2)
where $`R_{JQ}`$ is the distance from Jupiter to the quasar. The leading contribution to $`\mathrm{\Delta }(t_1,t_2)`$ is therefore

$$\mathrm{\Delta }(t_1,t_2)=\mathrm{\Delta }t_2\mathrm{\Delta }t_1=\frac{4G_NM_J}{c^3}\mathrm{ln}\left(\frac{\xi _1}{\xi _2}\right)=\frac{4G_NM_J\mathrm{\Delta }\xi }{\xi c^3}.$$
(3)
Let us determine $`\mathrm{\Delta }\xi =\xi _1\xi _2`$ in terms of $`\stackrel{}{B}`$. The electromagnetic rays that originate from the quasar are bent slightly as they pass by Jupiter by an amount $`\mathrm{\Delta }\phi `$ given by

$$\mathrm{\Delta }\phi =\frac{4G_NM_J}{\xi c^2}.$$

Eventually, one finds

$$\stackrel{}{n}\stackrel{}{B}=\mathrm{\Delta }\xi \left(1+\frac{4G_NM_JR_{EJ}}{\xi ^2c^2}\right)\mathrm{\Delta }\xi =\xi _1\xi _2,$$
(4)
because

$$\frac{4G_NM_JR_{EJ}}{\xi ^2c^2}\frac{4G_NM_JR_{EJ}}{R_J^2c^2}0.001,$$

where $`R_J`$ is the radius of Jupiter. In other words, within the solar system the angular deflection created by Jupiter can be neglected, and the separation between the rays remains essentially constant.

By substituting Eq.(4) into (3), one obtains the result for a static Jupiter

$$\mathrm{\Delta }(t_1,t_2)=\frac{4G_NM_J\stackrel{}{n}\stackrel{}{B}}{\xi c^3}.$$
(5)
When $`\xi =\theta _{obs}R_{EJ}`$ is used in Eq.(5), it reproduces the leading term in the notation of references A and B.

Let us now compute the $`v_J/c`$ corrections. This is simple to do by selecting an appropriate reference frame.

During the time in which the rays propagate from Jupiter to the Earth, Jupiter moves almost in a straight line with constant speed. In other words, the orbital motion of Jupiter around the Sun is not important. The same is true for the Earth. Therefore, observers on both planets can be considered as being inertial. Let us select an observational frame for which Jupiter is motionless. In this frame, the Earth appears to be moving with a velocity $`\stackrel{}{v}_E`$ equal to $`\stackrel{}{v}_J`$. Since Jupiter is not moving, Eq.(5) applies. However, the distance $`\stackrel{}{B}_{sf}`$ between points 1 and 2 as measured in the static frame is not equal to $`\stackrel{}{B}`$ as measured on Earth. Place a static observer at the point $`1`$ at time $`t_1`$ and another static observer at the point $`2`$ at time $`t_2`$. Have these observers make the time measurements. Then the situation is completely static and the formulas for the static case may be used.

During the time $`t_2t_1`$, the Earth moves a distance $`\stackrel{}{v}_E\left(t_2t_1\right)`$. Next, note that the leading contributions to $`t_2t_1`$ are

$$t_2t_1\frac{\stackrel{}{K}\stackrel{}{B}}{c}+\frac{\stackrel{}{n}\stackrel{}{B}\theta _{obs}}{c}+\mathrm{\Delta }(t_1,t_2),,$$
(6)
of which the first is the largest. Therefore,

$$\stackrel{}{B}_{sf}=\stackrel{}{B}\frac{\stackrel{}{K}\stackrel{}{B}}{c}\stackrel{}{v}_E+\frac{\stackrel{}{n}\stackrel{}{B}\theta _{obs}}{c}\stackrel{}{v}_E+\mathrm{\Delta }(t_1,t_2)\stackrel{}{v}_E.$$
(7)
The motion of Earth leads to two corrections to the static time delay difference in Eq.(5). Using $`\stackrel{}{n}\stackrel{}{B}_{sf}`$ in Eq.(5) leads to an additional term $`\frac{4G_NM_J\stackrel{}{n}\stackrel{}{v}_E\stackrel{}{K}\stackrel{}{B}}{\xi c^4}`$. The other correction arises if the Earth moves toward (or away from) Jupiter. In this case, the time delay is reduced (or increased) by the time $`\delta \mathrm{\Delta }(t_1,t_2)`$ it takes light to travel the distance determined by the difference between $`\stackrel{}{B}_{sf}`$ and $`\stackrel{}{B}`$. The corresponding correction due to the second and third terms of Eq.(7) is independent of $`G_N`$ and is a contribution to the first part of Eq.(1) that involves the difference in distances between the positions of the quasar and the observation points $`1`$ and $`2`$. The fourth term in Eq.(7) leads to

$$\delta \mathrm{\Delta }(t_1,t_2)=\frac{\stackrel{}{K}\stackrel{}{v}_E}{c}\mathrm{\Delta }(t_1,t_2)$$

One switches to the Earth frame using $`\stackrel{}{v}_E=\stackrel{}{v}_J`$. The final result is

$$\mathrm{\Delta }(t_1,t_2)=\frac{4G_NM_J}{\xi c^3}\left(\stackrel{}{n}\stackrel{}{B}\left(1+\frac{\stackrel{}{K}\stackrel{}{v}_J}{c}\right)+\frac{\stackrel{}{K}\stackrel{}{B}\stackrel{}{n}\stackrel{}{v}_J}{c}\right).$$
(8)
The correction factor $`\frac{\stackrel{}{K}\stackrel{}{v}_J}{c}`$ is present in references A and B. However, we find no $`1/\theta ^2`$ terms. In its place is the $`\stackrel{}{K}\stackrel{}{B}\stackrel{}{n}\stackrel{}{v}_J/c`$ term of Eq.(8).

Although the Shapiro time delay has effects created by the long-ranged gravitational force (e.g. see Eq.(2)), these effects cancel in the time difference of Eq.(1). In the static case, this is illustrated by Eq.(5), in which $`\mathrm{\Delta }(t_1,t_2)`$ is expressed in terms of the impact parameters of the electromagnetic waves, that is, quantities measurable in the vicinity of Jupiter. One therefore expects that long-ranged effects should not be present in $`\mathrm{\Delta }(t_1,t_2)`$ even in the non-static case. The $`1/\theta _{obs}^2`$ terms of reference A and B, however, grow with the Earth-Jupiter distance. On physical grounds, it seems unlikely that such terms are present, and our computation confirms this.

The leading term in $`\mathrm{\Delta }(t_1,t_2)`$ is of the order of $`100`$ nanoseconds for the Jupiter/quasar experiment. The $`v_J/c`$ corrections in Eq.(8) are at least 10,000 times smaller making them less than $`0.01`$ nanoseconds. It is therefore impossible that references  and  measured the speed of gravity. The $`v_J/c`$ corrections are also masked by larger corrections such as terms down by $`B/\xi `$ and $`\theta _{obs}`$, which are present but not shown in this Letter.

III. Comparison to References A and B

It is easy to find the source of the $`1/\theta ^2`$ effects in references A and B. In those works, the times $`s_1=t_1|\stackrel{}{x}_1\left(t_1\right)\stackrel{}{x}_J\left(s_1\right)|/c`$ and $`s_2=t_2|\stackrel{}{x}_2\left(t_2\right)\stackrel{}{x}_J\left(s_2\right)|/c`$ at which rays $`1`$ and $`2`$ pass by Jupiter are expanded in terms of the times $`t_1`$ and $`t_2`$ when the rays are observed at the points $`\stackrel{}{x}_1`$ and $`\stackrel{}{x}_2`$ on Earth. The differences between the $`s_i`$ and $`t_i`$ are sizeable, of order of $`R_{EJ}/c`$, and during this time, Jupiter moves a significant distance. See the dotted circle in Figure 1. References A and B define the angle $`\theta `$ in terms of the position of Jupiter at $`t_1`$. This is not the physically observed angle $`\theta _{obs}`$. For clarity, denote the angle of references A and B by $`\theta _{AB}`$.

When the electromagnetic waves from the quasar pass by Jupiter, sunlight that has been reflected off of Jupiter also heads toward Earth. Eventually, the various waves arrive on Earth. See Figure 1. It is evident that the angle $`\theta _{obs}`$ between the quasar and Jupiter observed by an astronomer on Earth is determined by Jupiter’s position at time $`s_1`$ and not $`t_1`$.

The reason for the $`1/\theta ^2`$ term in references A and B is due to the use of the artificial angle $`\theta _{AB}`$. The relation between $`\theta _{obs}`$ and $`\theta _{AB}`$ is

$$\theta _{obs}\theta _{AB}+\frac{\stackrel{}{n}\stackrel{}{v}_J}{c}.$$
(9)
When this result is substituted into the leading term of Eq.(5),

$$\mathrm{\Delta }(t_1,t_2)=\frac{4G_NM_J\stackrel{}{n}\stackrel{}{B}}{R_{EJ}\theta _{obs}c^3}=\frac{4G_NM_J\stackrel{}{n}\stackrel{}{B}}{R_{EJ}c^3}\left(\frac{1}{\theta _{AB}}\frac{\stackrel{}{n}\stackrel{}{v}_J}{c\theta _{AB}^2}\right),$$
(10)
and the $`1/\theta ^2`$ effect emerges.

References A and B express $`\mathrm{\Delta }(t_1,t_2)`$ as

$$\mathrm{\Delta }(t_1,t_2)=\frac{4G_NM_J}{c^3}\mathrm{ln}\left[\frac{r_{1J}\left(s_1\right)+\stackrel{}{K}\stackrel{}{r}_{1J}\left(s_1\right)}{r_{2J}\left(s_2\right)+\stackrel{}{K}\stackrel{}{r}_{2J}\left(s_2\right)}\right],$$
(11)
and then expands unwisely about $`t_1`$. The expansion is somewhat subtle since factors such as $`r_{1J}\left(s_1\right)+\stackrel{}{K}\stackrel{}{r}_{1J}\left(s_1\right)`$ are proportional to the small quantity $`\theta _1^2`$. A careful analysis reveals that Eqs.(8) and (29) of references A and B should have used

$$\mathrm{ln}\left[\frac{r_{1J}\left(s_1\right)+\stackrel{}{K}\stackrel{}{r}_{1J}\left(s_1\right)}{r_{2J}\left(s_2\right)+\stackrel{}{K}\stackrel{}{r}_{2J}\left(s_2\right)}\right]=\left[\frac{r_{1J}\left(t_1\right)+\stackrel{}{K}\stackrel{}{r}_{1J}\left(t_1\right)}{r_{2J}\left(t_1\right)+\stackrel{}{K}\stackrel{}{r}_{2J}\left(t_1\right)}\right]+\frac{2\stackrel{}{n}\stackrel{}{v}_J\stackrel{}{n}\stackrel{}{B}}{cr_{1J}\theta _{AB}^2}.$$
(12)
When Eq.(12) is substituted into Eq.(11), the $`1/\theta _{AB}^2`$ term of Eq.(10) due to the expansion in Eq.(9) is reproduced. This shows that our equations are consistent with the method used in references A and B.

Summarizing, (1) an analysis using a static frame allows one to easily compute the $`v_J/c`$ corrections from the static result and one finds no $`1/\theta _{obs}^2`$ terms, (2) physical considerations suggest that such terms are absent, (3) an ill-advised expansion about the arrival time $`t_1`$ in references A and B produces the artificial $`1/\theta _{AB}^2`$ effects, and finally, (4) consistency with references A and B is achieved after the expansion in Eq.(12) is used.

There is nothing wrong with using Eq.(10) as long as the instantaneous angle $`\theta _{AB}`$ of Figure 1 is used. However, in the analysis of data, reference  begins with Eq.(10) and replaces $`c`$ with $`c_g`$. This is a mistake since the difference between the angles $`\theta _{obs}`$ and $`\theta _{AB}`$ is due to the motion of Jupiter during the period in which light travels from Jupiter to the Earth and has nothing to do with gravity. It is not surprising therefore that when a fit to Eq.(10) is performed with $`c_g`$ replacing $`c`$ that the experimentally deduced value of $`c_g`$ is approximately $`c`$. Reference  is only measuring the leading contribution to the Shapiro time delay in rearranged form as must be the case since the $`v_J/c`$ effects were beyond detection with current radio telescopes.

IV. On the Notion of the Speed of Gravity

It is clear from the derivation using the static frame in Section II that the leading $`v_J/c`$ corrections involve the speed of light and not the speed of gravity, and there is a recent analysis that supports this claim. See also references  and . However, references A and B argue that $`v_J/c_g`$ should appear. The issue here is how does one extend Einstein’s general theory of relativity to allow the possibility that the speed of gravity $`c_g`$ is not equal to $`c`$. A reasonable approach is to assume that the effect of gravity propagates at $`c_g`$ instead of $`c`$. For example, in the retarded times and positions of Jupiter in formulas, one replaces $`v_J/c`$ by $`v_J/c_g`$. Hence, in the frame in which Jupiter is moving and the Earth is at rest, the $`v_J/c`$ effects are generated in the vicinity of Jupiter, and $`v_J/c_g`$ should appear in lieu of $`v_J/c`$ in Eq.(8) of $`\mathrm{\Delta }(t_1,t_2)`$. But consistancy demands that the computation of $`\mathrm{\Delta }(t_1,t_2)`$ be frame independent. Thus, there does not seem to be a consistant way to define the speed of gravity concept for the Jupiter/quasar experiment. In the static frame, the corrections are due to the speed of light, while in the Jupiter-moving frame they are due to the speed of gravity.

How then might one try to test $`c_gc`$ in Einstein’s theory of relativity? The static and Jupiter-moving frames are both inertial. If Jupiter happened to be accelerating toward (or away from) the quasar’s electromagnetic waves as they passed by the planet, then one would not be able to go back and forth between the two frames. Therefore, there is a reasonable chance that the speed of gravity concept could be defined for such a situation. The parameter $`c_g`$ would not be attached to velocity-dependent terms but to acceleration effects. Although it is worth exploring this possibility theoretically, it is unlikely that a system within or beyond our Solar System exists that generates an effect sufficiently large to be measurable with current instruments.

Acknowledgments

This work was supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, of the Department of Energy under contract number DE-AC03-76SF00098.

Figure Captions

Figure 1. The Motion of Electromagnetic Waves Relevant for the Jupiter/Quasar Experiment.
For clarity, the diagram is not drawn to scale.


format: text

text:
# References

SPIN FRUSTRATION IN AN EXACTLY SOLVABLE ISING-HEISENBERG DIAMOND CHAIN
Michal Jaščur and Jozef Strečka
Department of Theoretical Physics and Astrophysics, Faculty of Science,

P.J. Šafárik University, Moyzesova 16, 041 54 Košice, Slovak Republic

E-mail: jascur@kosice.upjs.sk, jozkos@pobox.sk

## Abstract

Exact results of a mixed spin-1/2 and spin-1 Ising-Heisenberg diamond chain are obtained applying a generalized decoration-iteration transformation. Depending on the values of external field, exchange parameters and anisotropy there exist six different phases in the ground state of the system. The magnetic order of these phases is briefly described in this work.
Keywords: Ising-Heisenberg diamond chain; Exact results

During last two decades low-dimensinal magnetic materials have become important subject in the research of magnetism. Among systems exhibiting interesting quantum magnetic phenomena one should mention those having a diamond-chain structure that have been recently intensively studied both theoretically and experimentally (see  and references therein). In order to describe these materials, the quantum anisotropic Heisenberg model has been frequently adopted. Though this model proves to be very useful, it brings mathematical complexities preventing one to obtain exact results even for 1D systems. For this reason, the main purpose of this work is to introduce a simplified version of the Heisenberg model which is exactly solvable for the diamond chain. We will investigate here the diamond chain consisting of Ising- (gray circles) and Heisenberg-type atoms (black circles), as it is depicted in Fig. 1. In order to obtain exact expression for the partition function of the model, we express the total Hamiltonian $`\widehat{}`$ as a sum of the bond Hamiltonians, i.e. $`\widehat{}=_{k=1}^N\widehat{}_k`$, where $`N`$ denotes the total number of the Ising-type atoms and the bond Hamiltonian $`\widehat{}_k`$ is given by (see Fig. 1.)

$`\widehat{}_k`$ $`=`$ $`J_1(\widehat{S}_{k1}^z+\widehat{S}_{k2}^z)(\widehat{\mu }_{k1}^z+\widehat{\mu }_{k2}^z)`$ (1)
$`+`$ $`J_2\left[\mathrm{\Delta }(\widehat{S}_{k1}^x\widehat{S}_{k2}^x+\widehat{S}_{k1}^y\widehat{S}_{k2}^y)+\widehat{S}_{k1}^z\widehat{S}_{k2}^z\right]`$
$``$ $`H_A(\widehat{\mu }_{k1}^z+\widehat{\mu }_{k2}^z)/2H_B(\widehat{S}_{k1}^z+\widehat{S}_{k2}^z)`$
In this equation, $`\widehat{S}_{k\alpha }^\gamma `$ and $`\widehat{\mu }_{k\alpha }^z`$ ($`\gamma =x,y,z;\alpha =1,2`$) represent relevant components of the standard spin-1 and spin-1/2 operators, respectively, the exchange interaction $`J_1`$ couples the nearest-neighbouring Ising and Heisenberg atoms, $`J_2`$ couples nearest-neighbouring Heisenberg pairs and $`\mathrm{\Delta }`$ represents the XXZ anisotropy parameter. Finally, the terms including $`H_A`$ nad $`H_B`$ describe interactions of relevant atoms with an external magnetic field. Next, introducing a generalized decoration-iteration transformation  and utilizing the commutability of the bond Hamiltonians, we easily find the following relation

$$𝒵(\beta ,J_1,J_2,\mathrm{\Delta },H_A,H_B)=A^N𝒵_0(\beta ,R,H_0),$$
(2)
where $`𝒵`$ represents the partition function of our model and $`𝒵_0`$ is the well-known partition function of a simple spin-1/2 Ising linear chain and $`\beta =1/(k_BT)`$. For unknown transformation parameters $`A`$, $`R`$ and $`H_0`$, we obtain relatively simple expressions :

$`A^4`$ $`=`$ $`(W_0+W_1)(W_0+W_2)(W_0+W_3)^2`$ (3)
$`\beta R`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{(W_0+W_1)(W_0+W_2)}{(W_0+W_3)^2}}`$ (4)
$`\beta H_0`$ $`=`$ $`\beta H_A\mathrm{ln}(W_0+W_1)+\mathrm{ln}(W_0+W_2)`$ (5)
$`W_0`$ $`=`$ $`\mathrm{exp}(t_2)+f(t_2/2,a,0,0)4`$ (6)
$`W_1`$ $`=`$ $`f(t_2,2t_1+2h_B,t_1+h_B,t_2\mathrm{\Delta })`$ (7)
$`W_2`$ $`=`$ $`f(t_2,2t_12h_B,t_1h_B,t_2\mathrm{\Delta })`$ (8)
$`W_3`$ $`=`$ $`f(t_2,2h_B,h_B,t_2\mathrm{\Delta })`$ (9)
with $`a=\sqrt{1+8\mathrm{\Delta }^2}`$, $`t_1=\beta J_1,t_2=\beta J_2,h_B=\beta H_B`$, and $`f(x_1,x_2,x_3,x_4)=2e^{x_1}\mathrm{cosh}x_2+4\mathrm{cosh}x_3\mathrm{cosh}x_4.`$ Eqs. (2)-(9) enable one to calculate all physical quantities of interest.

In order to illustrate some interesting properties of the system, we have depicted the ground-state phase diagram in the $`J_2H`$ space for the case of the antiferromagnetic isotropic model i.e., $`J_1>0`$, $`J_2>0`$, $`H_A=H_B=H`$ and $`\mathrm{\Delta }=1`$. To describe different phases, we have calculated and analysed the reduced magnetization of the Ising ($`m_i^z`$) and Heisenberg sublattice ( $`m_h^z`$) and various nearest-neighbour correlation functions ($`C_{hh}^{xx}`$, $`C_{hh}^{zz}`$, $`C_{ih}^{zz}`$ and $`C_{ii}^{zz}`$). All these quantities are defined in the standard way, for instance, $`m_i^z\widehat{\mu }_{k\alpha }^z,m_h^z\widehat{S}_{k\alpha }^z,C_{hh}^{xx}\widehat{S}_{k1}^x\widehat{S}_{k2}^x,`$ where $`\mathrm{}`$ means standard canonical average, subscripts $`i(h)`$ reffer to the Ising (Heisenberg) atoms and superscripts represent spatial components of the spin. As one can see from Fig.2, there exist six different phases separated by first-order phase transition boundaries. With the help of above introduced quantities we can characterize these phases as follows:
Ordered phase 1 (OP1):
$`m_i^z=1/2`$, $`m_h^z=1`$, $`C_{hh}^{zz}=1`$, $`C_{hh}^{xx}=0`$, $`C_{ih}^{zz}=1/2`$, $`C_{ii}^{zz}=1/4`$;
Ordered phase 2 (OP2):
$`m_i^z=1/2`$, $`m_h^z=1/2`$, $`C_{hh}^{zz}=0`$, $`C_{hh}^{xx}=1/2`$, $`C_{ih}^{zz}=1/4`$, $`C_{ii}^{zz}=1/4`$;
Ordered phase 3 (OP3):
$`m_i^z=1/2`$, $`m_h^z=0`$, $`C_{hh}^{zz}=2/3`$, $`C_{hh}^{xx}=2/3`$, $`C_{ih}^{zz}=0`$, $`C_{ii}^{zz}=1/4`$;
Ordered phase 4 (OP4):
$`m_i^z=1/2`$, $`m_h^z=1/2`$, $`C_{hh}^{zz}=0`$, $`C_{hh}^{xx}=1/2`$, $`C_{ih}^{zz}=1/4`$, $`C_{ii}^{zz}=1/4`$;
Saturated paramagnetic phase (SP):
$`m_i^z=1/2`$, $`m_h^z=1`$, $`C_{hh}^{zz}=1`$, $`C_{hh}^{xx}=0`$, $`C_{ih}^{zz}=1/4`$, $`C_{ii}^{zz}=1/4`$;
Frustrated phase (FP):
$`m_i^z=0`$, $`m_h^z=0`$, $`C_{hh}^{zz}=2/3`$, $`C_{hh}^{xx}=2/3`$, $`C_{ih}^{zz}=0`$, $`C_{ii}^{zz}=0`$.
As one can see the OP1 and SP represent standard phases that are usually observed also in the pure Ising models. On the other hand, the remaining four phases exhibit apparently a nontrivial spin ordering related to strong quantum fluctuations in the system. Due to this effect magnetic properties of the OP2 and OP4 can be explained within the valence-bond-solid picture as suggested in . Moreover, the stepwise magnetization curves exhibit plateaux at relevant ordered phases OP1, OP2, OP3 or OP4. In addition to the above mentioned phases, there appears for $`H/J=0`$ and $`J_2>J_1`$ an interesting frustrated phase (FP) having the nonzero entropy, namely, $`S/Nk_B=\mathrm{ln}2`$. Even without presenting further details here, it is clear that the system under consideration is very interesting and it will exhibit nontrivial properties also at finite temperatures.

Acknowledgment: This work was supported under grants APVT 20-009902 and VEGA 1/9034/02
FIGURE CAPTIONS:
Fig.1: Part of the of the mixed-spin diamond chain. Gray circles represent spin-1/2 Ising atoms and black ones spin-1 Heisenberg atoms. The ellipse demarcates $`k`$th bond (see text)
Fig.2: Ground-state phase diagram of the Ising-Heisenberg diamond chain in the $`J_2H`$ space. OP1, OP2, OP3, OP4, SP and FP denote different phases described in the text


text:
# Testing the statistical compatibility of independent data sets

## I Introduction

The essence of any scientific progress is the comparison of theoretical predictions to experimental data. Statistics provides the scientist with so-called goodness-of-fit tests, which allow to obtain well defined probability statements about the agreement of a theory with data. The by far most popular goodness-of-fit test dates back to 1900, when K. Pearson identified the minimum of a $`\chi ^2`$-function as a powerful tool to evaluate the quality of the fit pearson . However, it is known that the Pearson $`\chi _{\mathrm{min}}^2`$ test is not very restrictive in global analyses, where data from different experiments with a large number of data points are compared to a theory depending on many parameters. The reason for this is that in such a case a given parameter is often constrained only by a small subset of the data. If the rest of the data (which can contain many data points) are reasonably fitted, a possible problem in the fit of the given parameter is completely washed out by the large amount of data points. A discussion of this problem in various contexts can be found e.g. in Refs. collins ; Maltoni:2001bc ; solGof .

To evade this problem a modification of the original $`\chi _{\mathrm{min}}^2`$ test was proposed in Ref. Maltoni:2002xd  to evaluate the goodness-of-fit of neutrino oscillation data in the framework of four-neutrino models. There this method was called parameter goodness-of-fit (PG), and it can be applied when the global data consists of statistically independent subsets. The PG is based on parameter estimation and hence it avoids the problem of being diluted by many data points. It tests the compatibility of the different data sets in the framework of the given theoretical model. In this note we give a formal derivation of the probability distribution function (p.d.f.) for the test statistic of the PG, and discuss the application and interpretation of the PG on some examples. The original motivation for the PG was the analysis of neutrino oscillation data. However, the method may be very useful also in other fields of physics, especially where global fits of many parameters to data from several experiments are performed.

The outline of the paper is as follows. In Sec. II we define the PG and show that its construction is very similar to the one of the standard goodness-of-fit. The formal derivation of the p.d.f. for the PG test statistic is given in Sec. III, whereas in Sec. IV a discussion of the application and interpretation of the PG is presented. In Sec. V we consider the PG in the case of correlations due to theoretical errors, and we conclude in Sec. VI.

## II Goodness-of-fit tests

We would like to start the discussion by citing the goodness-of-fit definition given by the Particle Data Group (see Sec. 31.3.2. of Ref. pdg ): “Often one wants to quantify the level of agreement between the data and a hypothesis without explicit reference to alternative hypotheses. This can be done by defining a goodness-of-fit statistic, $`t`$, which is a function of the data whose value reflects in some way the level of agreement between the data and the hypothesis. \[…\] The hypothesis in question, say, $`H_0`$ will determine the p.d.f. $`g(t|H_0)`$ for the statistic. The goodness-of-fit is quantified by giving the $`p`$-value, defined as the probability to find $`t`$ in the region of equal or lesser compatibility with $`H_0`$ than the level of compatibility observed with the actual data. For example, if $`t`$ is defined such that large values correspond to poor agreement with the hypothesis, then the $`p`$-value would be

$$p=_{t_{\mathrm{obs}}}^{\mathrm{}}g(t|H_0)dt,$$
(1)
where $`t_{\mathrm{obs}}`$ is the value of the statistic obtained in the actual experiment.”

Let us stress that from this definition of goodness-of-fit one has complete freedom in choosing a test statistic $`t`$, as long as the correct p.d.f. for it is used.

### II.1 The standard goodness-of-fit

Consider $`N`$ random observables $`𝝂=(\nu _i)`$ and let $`\mu _i(𝜽)`$ denote the expectation value for the observable $`\nu _i`$, where $`𝜽=(\theta _\alpha )`$ are $`P`$ independent parameters which we wish to estimate from the data. Assuming that the covariance matrix $`S`$ is known one can construct the following $`\chi ^2`$-function:

$$\chi ^2(𝜽)=[𝝂𝝁(𝜽)]^TS^1[𝝂𝝁(𝜽)]$$
(2)
and use its minimum $`\chi _{\mathrm{min}}^2`$ as test statistic for goodness-of-fit evaluation:

$$t(𝝂)=\chi _{\mathrm{min}}^2.$$
(3)
The hypothesis we want to test determines the p.d.f. $`g(t)`$ for this statistics. Once the real experiments have been performed, giving the results $`𝝂_{\text{obs}}`$, the goodness-of-fit is given by the probability of obtaining a $`t`$ larger than $`t_{\mathrm{obs}}`$, as expressed by Eq. (1). We will refer to this procedure as standard goodness-of-fit (SG):

$$p_{\text{SG}}=_{\chi _{\mathrm{min}}^2(𝝂_{\mathrm{obs}})}^{\mathrm{}}g(t)dt.$$
(4)
The great success of this method is mostly due to a very powerful theorem, which was proven over 100 years ago by K. Pearson111Pearson uses the slightly different test statistic
$$\chi _{Pearson}^2=\underset{i}{}\frac{[\nu _i\mu _i(𝜽)]^2}{\mu _i(𝜽)}$$
and assumes that the $`\nu _i`$ are independent. We prefer to use instead the $`\chi ^2`$ of Eq. (2), because in this way also correlated data can be considered. pearson  and which greatly simplifies the task of calculating the integral in Eq. (4). It can be shown under quite general conditions (see e.g. Ref. cramer ) that $`\chi _{\mathrm{min}}^2`$ follows a $`\chi ^2`$-distribution with $`NP`$ degrees of freedom (d.o.f.), so that $`g(t)=f_{\chi ^2}(t,NP)`$. Therefore, the integral in Eq. (4) becomes:

$$p_{\text{SG}}=CL(\chi _{\mathrm{min}}^2(𝝂_{\mathrm{obs}}),NP)_{\chi _{\mathrm{min}}^2(𝝂_{\mathrm{obs}})}^{\mathrm{}}f_{\chi ^2}(t,NP)dt,$$
(5)
where $`CL(\chi ^2,n)`$ is the confidence level function (see e.g. Fig. 31.1 of Ref. pdg ).

In the following we propose a modification of the SG, for the case when the data can be divided into several statistically independent subsets.

### II.2 The parameter goodness-of-fit

Consider $`D`$ statistically independent sets of random observables $`𝝂^r=(\nu _i^r)`$ ($`r=1,\mathrm{},D`$), each consisting of $`N_r`$ observables ($`i=1,\mathrm{},N_r`$), with $`N_{\mathrm{tot}}=_rN_r`$. Now a theory depending on $`P`$ parameters $`𝜽=(\theta _\alpha )`$ is confronted with the data. The total $`\chi ^2`$ is given by

$$\chi _{\mathrm{tot}}^2(𝜽)=\underset{r=1}{\overset{D}{}}\chi _r^2(𝜽),$$
(6)
where

$$\chi _r^2(𝜽)=[𝝂^r𝝁^r(𝜽)]^TS_r^1[𝝂^r𝝁^r(𝜽)]$$
(7)
is the $`\chi ^2`$ of the data set $`r`$. Now we define

$$\overline{\chi }^2(𝜽)=\chi _{\mathrm{tot}}^2(𝜽)\underset{r=1}{\overset{D}{}}\chi _{r,\mathrm{min}}^2,$$
(8)
where $`\chi _{r,\mathrm{min}}^2=\chi _r^2(\widehat{}𝜽_r)`$, and $`\widehat{}𝜽_r(𝝂^r)`$ are the values of the parameters which minimize $`\chi _r^2`$. Instead of the total $`\chi ^2`$-minimum we propose now to use

$$t(𝝂)=\overline{\chi }_{\mathrm{min}}^2=\overline{\chi }^2(\stackrel{~}{}𝜽)$$
(9)
as test statistic for goodness-of-fit evaluation. In Eq. (9) $`\overline{\chi }_{\mathrm{min}}^2`$ is the minimum of $`\overline{\chi }^2`$ defined in Eq. (8), and $`\stackrel{~}{}𝜽`$ are the parameter values at the minimum of $`\overline{\chi }^2`$, or equivalently of $`\chi _{\mathrm{tot}}^2`$. If we now denote by $`\overline{g}(t)`$ the p.d.f. for this statistic, we can define the corresponding goodness-of-fit by means of Eq. (1), in complete analogy to the SG case:

$$p_{\text{PG}}=_{\overline{\chi }_{\mathrm{min}}^2(𝝂_{\mathrm{obs}})}^{\mathrm{}}\overline{g}(t)dt.$$
(10)
This procedure was proposed in Ref. Maltoni:2002xd  with the name parameter goodness-of-fit (PG). Its construction is very similar to the SG, except that now $`\overline{\chi }^2`$ rather than $`\chi ^2`$ is used to define the test statistic.

In the next section we will show that also in the case of the PG the calculation of the integral appearing in Eq. (10) can be greatly simplified. Let us define

$$P_r\text{rank}\left[\frac{𝝁^r}{𝜽}\right].$$
(11)
This corresponds to the number of independent parameters (or parameter combinations), constrained by a measurement of $`𝝁^r`$.222If in some pathological cases $`P_r`$ depends on the point in the parameter space Eq. (11) should be evaluated at the true values of the parameters, see Sec. III.2. Then under general condition $`\overline{\chi }_{\mathrm{min}}^2`$ is distributed as a $`\chi ^2`$ with $`P_c=_rP_rP`$ d.o.f., so that Eq. (10) reduces to:

$$p_{\text{PG}}=CL(\overline{\chi }_{\mathrm{min}}^2(𝝂_{\mathrm{obs}}),P_c).$$
(12)
## III The probability distribution function of $`\overline{\chi }_{\mathrm{min}}^2`$

In this section we derive the distribution of the test statistic for the PG. This can be done in complete analogy to the SG. Therefore, we start by reviewing the corresponding proof for the SG, see e.g. Ref. cramer .

### III.1 The standard goodness-of-fit

Let us start from the $`\chi ^2`$ defined in Eq. (2). Since the covariance matrix $`S`$ is a real, positive and symmetric matrix one can always find an orthogonal matrix $`O`$ and a diagonal matrix $`s`$ such that $`S^1=O^Ts^2O`$. Hence, we can write the $`\chi ^2`$ in the following way:

$$\chi ^2(𝜽)=[𝝂𝝁(𝜽)]^TS^1[𝝂𝝁(𝜽)]=𝐲(𝜽)^T𝐲(𝜽),$$
(13)
where we have defined the new variables $`𝐲(𝜽)=sO[𝝂𝝁(\theta )]`$. Let us denote the (unknown) true values of the parameters by $`𝜽^0`$ and we define

$$𝐱𝐲(𝜽^0)=sO[𝝂𝝁(𝜽^0)].$$
(14)
Now we assume that the $`x_i`$ are normal distributed with mean zero and the covariance matrix $`\mathrm{𝟏}_N`$, which in particular implies that they are statistically independent. This assumption is obviously correct if the data $`\nu _i`$ are normal distributed with mean $`\mu _i(𝜽^0)`$ and covariance matrix $`S`$. However, it can be shown (see e.g. Refs. cramer ; frodesen ; roe ) that this assumption holds for a large class of arbitrary p.d.f. for the data under quite general conditions, especially in the large sample limit, i.e. large $`\nu _i`$. Under this assumption it is evident that $`\chi ^2(𝜽^0)=𝐱^T𝐱`$ follows a $`\chi ^2`$-distribution with $`N`$ d.o.f.. According to Eq. (3) the test statistic $`t`$ for the SG is given by the minimum of Eq. (13). To derive the p.d.f. for $`t`$ we state the following proposition:

###### Proposition 1

Let $`\widehat{}𝛉`$ be the values of the parameters which minimize Eq. (13). Then

$$\chi _{\mathrm{min}}^2=\chi ^2(\theta ^0)\mathrm{\Delta }\chi ^2,$$
(15)
with $`\chi _{\mathrm{min}}^2=\widehat{𝐲}^T\widehat{𝐲}`$ and $`\widehat{𝐲}𝐲(\widehat{}𝛉)`$, has a $`\chi ^2`$-distribution with $`NP`$ d.o.f. and $`\mathrm{\Delta }\chi ^2`$ has a $`\chi ^2`$-distribution with $`P`$ d.o.f. and is statistically independent of $`\chi _{\mathrm{min}}^2`$.

A rigorous proof of this proposition is somewhat intricate and can be found e.g. in Ref. cramer . In the following we give an outline of the proof dispensing with mathematical details for the sake of clarity.

The $`\widehat{}𝜽`$ are obtained by solving the equations

$$\frac{\chi ^2}{\theta _\alpha }=2𝐲^T\frac{𝐲}{\theta _\alpha }=0.$$
(16)
It can be proved (see e.g. Ref. cramer ) under very general conditions that Eqs. (16) have a unique solution $`\widehat{}𝜽`$ which converges to the true values $`𝜽^0`$ in the large sample limit. In this sense it is a good approximation333Note that Eq. (17) is exact if the $`𝐲`$ depend linearly on the parameters $`𝜽`$. to write

$$\widehat{𝐲}𝐱+B(\widehat{}𝜽𝜽^0),$$
(17)
where we have defined the rectangular $`N\times P`$ matrix $`B`$ by

$$B\frac{𝐲}{𝜽}|_{𝜽^0}.$$
(18)
With out loss of generality we assume that444If $`\mathrm{rank}[B]=P^{}<P`$ some of the parameters $`\theta _\alpha `$ are not independent. In this case one can perform a change of variables and choose a new sets of parameters $`\theta _\beta ^{}`$, such that $`\chi ^2(𝜽^{})`$ depends only on the first $`P^{}`$ of them. The remaining parameters are not relevant for the problem and can be eliminated from the very beginning. When repeating the construction in the new set of variables, the number of parameters will be equal to the rank of $`B`$. $`\mathrm{rank}[B]=P`$. From Eq. (17) we obtain

$$\frac{𝐲}{𝜽}|_{\widehat{}𝜽}\frac{𝐲}{𝜽}|_{𝜽^0}=B.$$
(19)
Using this last relation in Eq. (16) we find that $`\widehat{𝐲}`$ fulfils $`\widehat{𝐲}^TB=0`$. Multiplying Eq. (17) from the left side by $`B^T`$ this leads to

$$B^T𝐱=B^TB(\widehat{}𝜽𝜽^0).$$
(20)
Using Eqs. (17) and (20) we obtain

$$\widehat{𝐲}^T\widehat{𝐲}=𝐱^T𝐱(\widehat{}𝜽𝜽^0)^TB^TB(\widehat{}𝜽𝜽^0).$$
(21)
The symmetric $`P\times P`$ matrix $`B^TB`$ can be written as $`B^TB=Rb^2R^T`$ with the orthogonal matrix $`R`$ and the diagonal matrix $`b`$, and Eq. (20) implies $`b^1R^TB^T𝐱=bR^T(\widehat{}𝜽𝜽^0)`$. Defining the $`N\times P`$ matrix

$$HBRb^1$$
(22)
we find $`(\widehat{}𝜽𝜽^0)^TB^TB(\widehat{}𝜽𝜽^0)=𝐱^THH^T𝐱`$, and Eq. (21) becomes

$$\widehat{𝐲}^T\widehat{𝐲}=𝐱^T(\mathrm{𝟏}_NHH^T)𝐱.$$
(23)
Note that the matrix $`H`$ obeys the orthogonality relation $`H^TH=\mathrm{𝟏}_P`$, showing that the $`P`$ column vectors of length $`N`$ in $`H`$ are orthogonal. We can add $`NP`$ columns to the matrix $`H`$ completing it to an orthogonal $`N\times N`$ matrix: $`V=(H,K)`$. Here $`K`$ is an $`N\times (NP)`$ matrix with $`K^TK=\mathrm{𝟏}_{(NP)}`$, $`H^TK=0`$ and the completeness relation

$$VV^T=HH^T+KK^T=\mathrm{𝟏}_N.$$
(24)
Now we transform to the new variables

$$𝐱^{}=V^T𝐱,𝐱^{}=\left(\begin{array}{c}𝐯\\ 𝐰\end{array}\right)=\left(\begin{array}{c}H^T𝐱\\ K^T𝐱\end{array}\right),$$
(25)
where $`𝐯=H^T𝐱`$ is a vector of length $`P`$ and $`𝐰=K^T𝐱`$ is a vector of length $`NP`$. In general, if the covariance matrix of the random variables $`𝐱`$ is $`S`$, then the covariance matrix $`S^{}`$ of $`𝐱^{}=V^T𝐱`$ is given by $`S^{}=V^TSV`$. Hence, since in the present case the $`x_i`$ are normal distributed with mean 0 and covariance matrix $`\mathrm{𝟏}_N`$ the same is true for the $`x_i^{}`$. In particular also $`𝐯`$ and $`𝐰`$ are statistically independent. Using Eqs. (23) and (24) we deduce

$$\widehat{𝐲}^T\widehat{𝐲}=𝐱^T(\mathrm{𝟏}_NHH^T)𝐱=𝐱^TKK^T𝐱=𝐰^T𝐰$$
(26)
proving that $`\chi _{\mathrm{min}}^2=\widehat{𝐲}^T\widehat{𝐲}`$ has a $`\chi ^2`$-distribution with $`NP`$ d.o.f.. Finally, we obtain

$$\mathrm{\Delta }\chi ^2=\chi ^2(\theta ^0)\chi _{\mathrm{min}}^2=𝐱^T𝐱\widehat{𝐲}^T\widehat{𝐲}=𝐱^THH^T𝐱=𝐯^T𝐯,$$
(27)
showing that $`\mathrm{\Delta }\chi ^2`$ has a $`\chi ^2`$-distribution with $`P`$ d.o.f. and is statistically independent of $`\chi _{\mathrm{min}}^2`$. $`\mathrm{}`$

### III.2 The parameter goodness-of-fit

Moving now to the PG we generalize in an obvious way the formalism of the previous section by attaching and index $`r`$ for the data set to each quantity. We have

$$\chi _{\mathrm{tot}}^2(𝜽)=\underset{r}{}𝐲_r^T(𝜽)𝐲_r(𝜽),\chi _{\mathrm{tot}}^2(𝜽^0)=\underset{r}{}𝐱_r^T𝐱_r,$$
(28)
and

$$\overline{\chi }^2(𝜽)\chi _{\mathrm{tot}}^2(𝜽)\underset{r}{}\chi _{\mathrm{r},\mathrm{min}}^2=\underset{r}{}\left[𝐲_r^T(𝜽)𝐲_r(𝜽)\widehat{𝐲}_r^T\widehat{𝐲}_r\right].$$
(29)
###### Proposition 2

Let $`\stackrel{~}{}𝛉`$ be the values of the parameters which minimize $`\overline{\chi }^2(𝛉)`$, or equivalently $`\chi _{\mathrm{tot}}^2(𝛉)`$. Then $`\overline{\chi }_{\mathrm{min}}^2=\overline{\chi }^2(\stackrel{~}{}𝛉)`$ follows a $`\chi ^2`$-distribution with $`P_c`$ d.o.f., with

$$P_c𝒫P,𝒫\underset{r=1}{\overset{D}{}}P_r,P_r\mathrm{rank}[B_r]\text{and}B_r\frac{𝐲_r}{𝜽}|_{𝜽^0}.$$
(30)
The matrices $`B_r`$ are of order $`N_r\times P`$. Since a given data set $`r`$ may depend only on some of the $`P`$ parameters, or on some combination of them, in general one has to consider the possibility of $`P_rP`$.555Note that the definition of $`P_r`$ in Eq. (30) is equivalent to the one given in Eq. (11). This means that the symmetric $`P\times P`$ matrix $`B_r^TB_r`$ can be writen as $`R_rb_r^Tb_rR_r^T`$, where $`R_r`$ is an orthogonal matrix and $`b_r`$ is a $`P_r\times P`$ “diagonal” matrix, such that the diagonal $`P\times P`$ matrix $`b_r^Tb_r`$ will have $`P_r`$ non-zero entries. Let us now define the $`P\times P_r`$ “diagonal” matrix $`b_r^1`$ in such a way that $`(b_r^1)_{ii}1/(b_r)_{ii}`$ for each of the $`P_r`$ non-vanishing entries of $`b_r`$, and all other elements are zero. In analogy to Eq. (22) we introduce now the matrices

$$H_rB_rR_rb_r^1,$$
(31)
which are of order $`N_r\times P_r`$. To prove Proposition 2 we define the vectors of length $`𝒫`$

$$𝐘(𝜽)\left(\begin{array}{c}H_1^T𝐲_1(𝜽)\\ \mathrm{}\\ H_D^T𝐲_D(𝜽)\end{array}\right),𝐗\left(\begin{array}{c}H_1^T𝐱_1\\ \mathrm{}\\ H_D^T𝐱_D\end{array}\right)=\left(\begin{array}{c}𝐯_1\\ \mathrm{}\\ 𝐯_D\end{array}\right).$$
(32)
In the first part of the proof we show that $`\overline{\chi }_{\mathrm{min}}^2=\stackrel{~}{𝐘}^T\stackrel{~}{𝐘}`$ with $`\stackrel{~}{𝐘}𝐘(\stackrel{~}{}𝜽)`$. With arguments similar to the ones leading to Eq. (21) we find

$$\underset{r}{}\stackrel{~}{𝐲}_r^T\stackrel{~}{𝐲}_r=\underset{r}{}𝐱_r^T𝐱_r(\stackrel{~}{}𝜽𝜽^0)^T\underset{r}{}B_r^TB_r(\stackrel{~}{}𝜽𝜽^0).$$
(33)
Using further Eq. (23) for each $`r`$ we obtain

$`\overline{\chi }_{\mathrm{min}}^2`$ $`=`$ $`{\displaystyle \underset{r}{}}\stackrel{~}{𝐲}_r^T\stackrel{~}{𝐲}_r{\displaystyle \underset{r}{}}\widehat{𝐲}_r^T\widehat{𝐲}_r`$ (34)
$`=`$ $`{\displaystyle \underset{r}{}}𝐱_r^TH_rH_r^T𝐱_r(\stackrel{~}{}𝜽𝜽^0)^T{\displaystyle \underset{r}{}}B_r^TB_r(\stackrel{~}{}𝜽𝜽^0).`$
On the other hand we can use that the minimum values $`\stackrel{~}{}𝜽`$ are converging to the true values $`𝜽^0`$ in the large sample limit and write $`\stackrel{~}{𝐘}𝐗+(\stackrel{~}{}𝜽𝜽^0)`$, where we have defined the $`𝒫\times P`$ matrix

$$\frac{𝐘}{𝜽}|_{𝜽^0}=\left(\begin{array}{c}H_1^TB_1\\ \mathrm{}\\ H_D^TB_D\end{array}\right).$$
(35)
Without loss of generality we assume that $`\mathrm{rank}[]=P`$. Again, with arguments similar to the ones leading to Eq. (21) we derive

$$\stackrel{~}{𝐘}^T\stackrel{~}{𝐘}=𝐗^T𝐗(\stackrel{~}{}𝜽𝜽^0)^T^T(\stackrel{~}{}𝜽𝜽^0).$$
(36)
Using Eq. (31) it is easy to show that $`^T=_rB_r^TB_r`$, and by comparing Eqs. (36) and (34) we can readily verify the relation $`\overline{\chi }_{\mathrm{min}}^2=\stackrel{~}{𝐘}^T\stackrel{~}{𝐘}`$.

To complete the proof we identify $`𝐘𝐲`$ and $`𝐗𝐱`$ and proceed in perfect analogy to the proof of Proposition 1 given in Sec. III.1. In particular, from the arguments presented there it follows that the elements of $`𝐯_r`$ are $`P_r`$ independent Gaussian variables with mean zero and variance one. Since the $`D`$ data sets are assumed to be statistically independent the vector $`𝐗`$ contains $`𝒫`$ independent Gaussian variables with mean zero and variance one. In analogy to the matrices $`H,K`$ of Sec. III.1 we obtain now the $`𝒫\times P`$ matrix $``$ and the $`𝒫\times P_c`$ matrix $`𝒦`$, which fulfil $`^T+𝒦𝒦^T=\mathrm{𝟏}_𝒫`$, and Eq. (36) becomes

$$\stackrel{~}{𝐘}^T\stackrel{~}{𝐘}=𝐗^T(\mathrm{𝟏}_𝒫^T)𝐗.$$
(37)
In analogy to the vector $`𝐰`$ from Eq. (25) we define now $`𝐖𝒦^T𝐗`$, containing $`P_c=𝒫P`$ independent Gaussian variables with mean zero and variance one, and Eq. (37) gives

$$\stackrel{~}{𝐘}^T\stackrel{~}{𝐘}=𝐗^T𝒦𝒦^T𝐗=𝐖^T𝐖.$$
(38)
From Eq. (38) it is evident that $`\overline{\chi }_{\mathrm{min}}^2=\stackrel{~}{𝐘}^T\stackrel{~}{𝐘}`$ follows a $`\chi ^2`$-distribution with $`P_c`$ d.o.f.. $`\mathrm{}`$

Let us conclude this section by noting that both Proposition 1 and 2 are exact if the data are multi-normally distributed and the theoretical predictions $`𝝁`$, $`𝝁^r`$ depend linearly on the parameters $`𝜽`$. If these requirements are not fulfilled the simplified expressions (5) and (12) are valid only approximately, and to calculate the SG and the PG one should in principle use the general formulas (4) and (10) instead. However, we want to stress that under rather general conditions $`\chi _{\mathrm{min}}^2`$ and $`\overline{\chi }_{\mathrm{min}}^2`$ will be distributed as a $`\chi ^2`$ in the large sample limit (i.e. for large $`𝝂`$ and $`𝝂^r`$, respectively), so that even in the general case Eqs. (5) and (12) can still be used.

## IV Examples and discussion

In this section we illustrate the application of the PG on some examples. In Sec. IV.1 we show that in the simple case of two measurements of a single parameter the PG is identical to the intuitive method of considering the difference of the two measurements, and in Sec. IV.2 we show the consistency of the PG and the SG in the case of independent data points. In Sec. IV.3 we discuss the application of the PG to neutrino oscillation data in the framework of a sterile neutrino scheme. This problem was the original motivation to introduce the PG in Ref. Maltoni:2002xd . In Sec. IV.4 we add some general remarks on the PG.

### IV.1 The determination of one parameter by two experiments

Let us consider two data sets observing the data points $`𝝂^1=(\nu _i^1)`$ ($`i=1,\mathrm{},N_1`$) and $`𝝂^2=(\nu _i^2)`$ ($`i=1,\mathrm{},N_2`$). Further, we assume that the expectation values for both data sets can be calculated from a theory depending on one parameter $`\eta `$: $`𝝁^r(\eta )`$ ($`r=1,2`$), and all $`\nu _i^r`$ are independent and normal distributed around the expectation values with variance $`\sigma _i^r`$. Then we have the following $`\chi ^2`$-functions for the two data sets $`r=1,2`$:

$$\chi _r^2(\eta )=\underset{i=1}{\overset{N_r}{}}\left(\frac{\nu _i^r\mu _i^r(\eta )}{\sigma _i^r}\right)^2=\chi _{r,\mathrm{min}}^2+\left(\frac{\widehat{\eta }_r\eta }{\widehat{\sigma }_r}\right)^2,$$
(39)
where $`\widehat{\eta }_r=\widehat{\eta }_r(𝝂^r)`$ is the value of the parameter at the $`\chi ^2`$-minimum of data set $`r`$. Now one may ask the question whether the results of the two experiments are consistent. More precisely, we are interested in the probability to obtain $`\widehat{\eta }_1`$ and $`\widehat{\eta }_2`$ under the assumption that both result from the same true value $`\eta ^0`$.

A standard method (see e.g. Ref. frodesen  Sec. 14.3) to answer this question is to consider the variable

$$z=\frac{\widehat{\eta }_1\widehat{\eta }_2}{\sqrt{\widehat{\sigma }_1^2+\widehat{\sigma }_2^2}}.$$
(40)
If the theory is correct $`z`$ is normal distributed with mean zero and variance one. Hence we can answer the question raised above by citing the probability to obtain $`|z||z_{\mathrm{obs}}|`$:

$$p=1_{|z_{\mathrm{obs}}|}^{|z_{\mathrm{obs}}|}f_N(z;0,1)dz,$$
(41)
where $`f_N`$ denotes the normal distribution.

If the PG is applied to this problem, one obtains from Eq. (39)

$$\overline{\chi }^2(\eta )=\left(\frac{\widehat{\eta }_1\eta }{\widehat{\sigma }_1}\right)^2+\left(\frac{\widehat{\eta }_2\eta }{\widehat{\sigma }_2}\right)^2,$$
(42)
and after some simple algebra one finds $`\overline{\chi }_{\mathrm{min}}^2=z^2`$, where $`z`$ is given in Eq. (40). Obviously, applying Eq. (12) to calculated the $`p`$-value according to the PG with the relevant number of d.o.f. $`P_c=21=1`$ leads to the same result as Eq. (41).

Hence, we arrive at the conclusion that in this simple case of testing the compatibility of two measurements for the mean of a Gaussian, the PG is identical to the intuitive method of testing whether the difference of the two values is consistent with zero.

### IV.2 Consistency of PG and SG for independent data points

As a further example of the consistency of the PG method we consider the case of $`N`$ statistically independent data points $`\nu _i`$. Let us denote by $`\sigma _i`$ the standard deviation of the observation $`\nu _i`$ ($`i=1,\mathrm{},N`$), and the corresponding theoretical prediction by $`\mu _i(𝜽)`$, where $`𝜽`$ is the vector of $`P`$ parameters. For simplicity, we assume that each of the $`\mu _i`$ depends at least on one parameter. Then the $`\chi ^2`$ is given by

$$\chi ^2(𝜽)=\underset{i=1}{\overset{N}{}}\chi _i^2(𝜽),\text{where}\chi _i^2(𝜽)=\frac{[\nu _i\mu _i(𝜽)]^2}{\sigma _i^2},$$
(43)
and from the SG construction (see Sec. III.1) we know that $`\chi _{\mathrm{min}}^2`$ follows a $`\chi ^2`$-distribution with $`NP`$ degrees of freedom. On the other hand, if we consider each single data point as an independent data set and we apply the PG construction, we easily see that $`\chi _{i,\mathrm{min}}^2=0`$ for each $`i`$. This implies $`\overline{\chi }^2(𝜽)=\chi ^2(𝜽)`$, and in particular $`\overline{\chi }_{\text{min}}^2=\chi _{\text{min}}^2`$. Therefore, for the specific case considered here one expects that SG and PG are identical.

To show that this is really the case let us first note that each matrix $`\mu _i/𝜽`$ consists just of *a single line*, and therefore it obviously has rank one. Hence, Eq. (11) gives $`P_i=1`$ for each $`i`$. This reflects the fact that from the measurement of a single observable we cannot derive *independent* bounds on $`P`$ parameters, but only *a single combination* of them is constrained. Therefore, the number of d.o.f. relevant for the calculation of the PG is given by $`P_c=_{i=1}^NP_iP=NP`$, which is exactly the number of d.o.f. relevant for the SG. Hence, we have shown that in the considered case the two methods are equivalent and consistent.

### IV.3 Application to neutrino oscillation data

In this section we use real data from neutrino oscillation experiments to discuss the application of the PG and to compare it to the SG. We consider the so-called (2+2) neutrino mass scheme, where a fourth (sterile) neutrino is introduced in addition to the three standard model neutrinos. In general this model is characterized by 9 parameters: 3 neutrino mass-squared differences $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$, $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$, $`\mathrm{\Delta }m_{\mathrm{LSND}}^2`$ and 6 mixing parameters $`\theta _{\mathrm{sol}}`$, $`\theta _{\mathrm{atm}}`$, $`\theta _{\mathrm{LSND}}`$, $`d_\mu `$, $`\eta _s`$, $`\eta _e`$. The interested reader can find precise definitions of the parameters, applied approximations, an extensive discussion of physics aspects, and references in Refs. Maltoni:2002xd ; Maltoni:2001bc ; Maltoni:2002ni . Here we are interested mainly in the statistical aspects of the analysis, and therefore we consider a simplified scenario.

We do not include LSND, KARMEN and all the experiments sensitive to $`\mathrm{\Delta }m_{\mathrm{LSND}}^2`$ and the corresponding mixing angle $`\theta _{\mathrm{LSND}}`$. Hence, we are left with three data sets from solar, atmospheric and reactor neutrino experiments. The solar data set includes the current global solar neutrino data from the SNO, Super-Kamiokande, Gallium and Chlorine experiments, making a total of 81 data points, whereas the atmospheric data sample includes 65 data points from the Super-Kamiokande and MACRO experiments (for details of the solar and atmospheric analysis see Ref. Maltoni:2002ni ). In the reactor data set we include only the data from the KamLAND and the CHOOZ experiments, leading to a total of $`13+14=27`$ data points Maltoni:2002aw ; Apollonio:2002gd . In general the reactor experiments (especially CHOOZ) depend in addition to $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ and $`\theta _{\mathrm{sol}}`$ also on $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ and a further mixing parameter $`\eta _e`$. However, we adopt here the approximation $`\eta _e=1`$, which is very well justified in the (2+2) scheme Maltoni:2001bc . This implies that the dependence on $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ disappears and we are left with the parameters $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ and $`\theta _{\mathrm{sol}}`$ for both reactor experiments, KamLAND as well as CHOOZ.

Under these approximations the experimental data sets we are using are described only by the 6 parameters $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$, $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$, $`\theta _{\mathrm{sol}}`$, $`\theta _{\mathrm{atm}}`$, $`\eta _s`$, $`d_\mu `$. The parameter structure is illustrated in Fig. 1. This simplified analysis serves well for discussing the statistical aspects of the problem; a more general treatment including a detailed discussion of the physics is given in Refs. Maltoni:2001bc ; Maltoni:2002xd . In Tab. 1 we summarize the parameter dependence, the number of data points, the minimum values of the $`\chi ^2`$-functions and the resulting SG. We observe that all the data sets analyzed alone give a very good fit. Let us remark that especially in the case of reactor and atmospheric data the SG is suspicious high. This may indicate that the errors have been estimated very conservatively.

In Tab. 2 we show the results of an SG and PG analysis for various combinations of the three data sets. In the first three lines in the table only two out of the three data sets are combined. By combining solar and atmospheric neutrino data we find a $`\chi _{\mathrm{min}}^2`$ of 126.7. With the quite large number of d.o.f. of 140 this gives an excellent SG of 78.3%. If however, the PG is applied we obtain a goodness-of-fit of only $`3.54\times 10^6`$. The reason for this very bad fit can be understood from Figs. 1 and 2. From Fig. 1 one finds that solar and atmospheric data are coupled by the parameter $`\eta _s`$. In Fig. 2 the $`\mathrm{\Delta }\chi ^2`$ is shown for both sets as a function of this parameter. We find that there is indeed significant disagreement between the two data sets:666The physical reason for this is that both data sets strongly disfavour oscillations into sterile neutrinos. Since it is a generic prediction of the (2+2) scheme that the sterile neutrino must show up either in solar or in atmospheric neutrino oscillations the model is ruled out by the PG test Maltoni:2002xd . solar data prefers values of $`\eta _s`$ close to zero, whereas atmospheric data prefers values close to one. There are two reasons why this strong disagreement does not show up in the SG. First, since the SG of both data sets alone is very good, there is much room to “hide” some problems in the combined analysis. Second, because of the large number of data points many of them actually might not be sensitive to the parameter $`\eta _s`$, where the disagreement becomes manifest. Hence, the problem in the combined fit becomes diluted due to the large number of data points. We conclude that the PG is very sensitive to disagreement of the data sets, even in cases where the individual $`\chi ^2`$-minima are very low, and when the number of data points is large.

In the reactor + solar analysis one finds complete agreement between the two data sets for the SG as well as for the PG. This reflects the fact that the determination of the parameters $`\theta _{\mathrm{sol}}`$ and $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ from reactor and solar neutrino experiment are in excellent agreement Maltoni:2002aw . Finally, in the case of the combined analysis of reactor and atmospheric data the PG cannot be applied. In our approximation these data sets have no parameter in common as one can see in Fig. 1. Hence, it makes no sense to test their compatibility, or even to combine them at all.

In the lower part of Tab. 2 we show the results from combining all three data sets. By comparing these results with the one from the solar + atmospheric analysis one can appreciate the advantage of the PG. If we add only the 13 data points from KamLAND to the solar and atmospheric samples we observe that the SG improves from 78.3% to 88.1%, whereas if both reactor experiments are included we obtain an SG of 95.0%. This demonstrates that the SG strongly depends on the number of data points. Especially the 14 data points from CHOOZ contain nearly no relevant information, since the best fit values of $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ and $`\theta _{\mathrm{sol}}`$ are in the no-oscillation regime for CHOOZ implying that the $`\chi ^2`$ is flat. Moreover, since reactor data are not sensitive to the parameter $`\eta _s`$ (see Fig. 1) the disagreement between solar and atmospheric data becomes even more diluted by the additional reactor data points. This clearly illustrates that the SG can be drastically improved by adding data which contains no information on the relevant parameters. Also the PG improves slightly by adding reactor data, reflecting the good agreement between solar and reactor data. However, the resulting PG still is very small due to the disagreement between solar and atmospheric data in the model under consideration. Moreover, the PG is completely unaffected by the addition of the CHOOZ data, because the $`\chi ^2`$ of CHOOZ is flat in the relevant parameter region, and the PG is sensitive only to the parameter dependence of the data sets.

Finally, we mention that in view of the analyses shown in Tab. 2 the meaning of $`P_c`$, the number of d.o.f. for PG becomes clear. It corresponds to the number of parameters coupling the data sets. Solar and atmospheric data are coupled only by $`\eta _s`$, hence $`P_c=1`$, whereas reactor and solar data are coupled by $`\theta _{\mathrm{sol}}`$ and $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ and $`P_c=2`$. Atmospheric data has no parameter in common with reactor data, therefore $`P_c=0`$. In the combination of reactor + solar + atmospheric data sets the three parameters $`\eta _s,\theta _{\mathrm{sol}},\mathrm{\Delta }m_{\mathrm{sol}}^2`$ provide the coupling and $`P_c=3`$.

### IV.4 General remarks on the PG

(a) Using the relation $`\overline{\chi }_{\mathrm{min}}^2=_r\mathrm{\Delta }\chi _r^2(\stackrel{~}{}𝜽)`$ one can obtain more insight into the quality of the fit by considering the contribution of each data set to $`\overline{\chi }_{\mathrm{min}}^2`$. If the PG is poor it is possible to identify the data sets leading to the problems in the fit by looking at the individual values of $`\mathrm{\Delta }\chi _r^2(\stackrel{~}{}𝜽)`$. In this sense the PG is similar to the so-called “pull approach” discussed in Ref. Fogli:2002pt  in relation with solar neutrino analysis.

(b) One should keep in mind that the PG is completely insensitive to the goodness-of-fit of the individual data sets. Because of the subtraction of the $`\chi _{r,\mathrm{min}}^2`$ in Eq. (8) all the information on the quality of the fit of the data sets alone is lost. One may benefit from this property if the SG of the individual data sets is very good (see the example in Sec. IV.3). On the other hand, if e.g. one data set gives a bad fit on its own this will not show up in the PG. Only the compatibility of the data sets is tested, irrespective of their individual SGs.

(c) The PG might be also useful if one is interested whether a data set consisting of very few data points is in agreement with a large data sample.777For example one could think of a combination of the 19 neutrino events from the Super Nova 1987A with the high statistics global solar neutrino data. The SG of the combined analysis will be completely dominated by the large sample and the information contained in the small data sample may be drowned out by the large number of data points. In such a case the PG can give valuable information on the compatibility of the two sets, because it is not diluted by the number of data points in each set and it is sensitive only to the parameter dependence of the sets.

## V Correlations due to theoretical uncertainties

One of the limitations of the PG is that it can be applied only if the data sets are statistically independent. In many physically interesting situations (for example, different solar neutrino experiments) this is not the case since theoretical uncertainties introduce correlations between the results of different – and otherwise independent – experiments. However, in such a case one can take advantage of the so-called pull approach, which, as demonstrated in Ref. Fogli:2002pt , is equivalent to the usual covariance method. In that paper it was shown that if correlations due to theoretical errors exist, it is possible to account for them by introducing new parameters $`\xi _a`$ and adding penalty functions to the $`\chi ^2`$. In this way it is possible to get rid of unwanted correlations and the PG can be applied. The correlation parameters $`\xi _a`$ should be treated in the same way as the parameters $`𝜽`$ of the theoretical model.

In this section we illustrate this procedure by considering a generic experiment with an uncertainty on the normalisation of the predicted number of events. Let the experiment observe some energy spectrum which is divided into $`N`$ bins. The theoretical prediction for the bin $`i`$ is denoted by $`\mu _i(𝜽)`$ depending on $`P`$ parameters $`𝜽`$. In praxis often $`\mu _i(𝜽)`$ is not known exactly. Let us consider the case of a fully correlated relative error $`\sigma _{\mathrm{th}}`$. A common method to treat such an error is to add statistical and theoretical errors in quadrature, leading to the correlation matrix

$$S_{ij}(𝜽)=\delta _{ij}\sigma _{i,\mathrm{stat}}^2+\sigma _{\mathrm{th}}^2\mu _i(𝜽)\mu _j(𝜽),$$
(44)
where $`\sigma _{i,\mathrm{stat}}`$ is the statistical error in the bin $`i`$. In the case of neutrino oscillation experiments such a correlated error results e.g. from an uncertainty of the initial flux normalisation or of the fiducial detector volume. The $`\chi ^2`$ is given by

$$\chi ^2(𝜽)=\underset{i,j=1}{\overset{N}{}}[\nu _i\mu _i(𝜽)]S_{ij}^1(𝜽)[\nu _j\mu _j(𝜽)],$$
(45)
where $`\nu _i`$ are the observations. As shown in Ref. Fogli:2002pt , instead of Eq. (45) we can equivalently use

$$\chi ^2(𝜽,\xi )=\underset{i=1}{\overset{N}{}}\left(\frac{\nu _i\xi \mu _i(𝜽)}{\sigma _{i,\mathrm{stat}}}\right)^2+\left(\frac{\xi 1}{\sigma _{\mathrm{th}}}\right)^2,$$
(46)
and minimize with respect to the new parameter $`\xi `$.

On the other hand, if $`\xi `$ is considered as an additional parameter, on the same footing as $`𝜽`$, all the data points are formally uncorrelated and it is straight forward to apply the PG. Subtracting the minimum of the first term in Eq. (46) with respect to $`𝜽`$ and $`\xi `$ one obtains

$$\overline{\chi }^2(𝜽,\xi )=\mathrm{\Delta }\chi ^2(𝜽,\xi )+\left(\frac{\xi 1}{\sigma _{\mathrm{th}}}\right)^2.$$
(47)
The external information on the parameter $`\xi `$ represented by the second term in Eq. (47) is considered as an additional data set. Evaluating the minimum of Eq. (47) for 1 d.o.f. is a convenient method to test if the best fit point of the model is in agreement with the constraint on the over-all normalisation. In particular one can identify whether a problem in the fit comes from the spectral shape (first term) or the total rate (second term).

Moreover, one may like to divide the data into two parts, set I consisting of bins $`1,\mathrm{},n`$ and and set II consisting of bins $`n,\mathrm{},N`$, and test whether these data sets are compatible. Eq. (46) can be written as

$$\chi ^2(𝜽,\xi )=\underset{i=1}{\overset{n}{}}\left(\frac{\nu _i\xi \mu _i(𝜽)}{\sigma _{i,\mathrm{stat}}}\right)^2+\underset{i=n}{\overset{N}{}}\left(\frac{\nu _i\xi \mu _i(𝜽)}{\sigma _{i,\mathrm{stat}}}\right)^2+\left(\frac{\xi 1}{\sigma _{\mathrm{th}}}\right)^2,$$
(48)
and subtracting the minima of the two first terms gives the $`\overline{\chi }^2`$ relevant for the PG:

$$\overline{\chi }^2(𝜽,\xi )=\mathrm{\Delta }\chi _\mathrm{I}^2(𝜽,\xi )+\mathrm{\Delta }\chi _{\mathrm{II}}^2(𝜽,\xi )+\left(\frac{\xi 1}{\sigma _{\mathrm{th}}}\right)^2.$$
(49)
Assuming that the data sets I and II both depend on all $`P`$ parameters $`𝜽`$ the minimum of this $`\overline{\chi }^2`$ has to be evaluated for $`P+2`$ d.o.f. to obtain the PG. This procedure tests whether the data sets I and II are consistent with each other and the constraint on the over-all normalisation. By considering the relative contributions of the three terms in Eq. (49) it is possible to identify potential problems in the fit. For example one may test whether a bad fit is dominated only by a small subset of the data, e.g. a few bins at the low or high end of the spectrum. Alternatively, the two data sets I and II can come from two different experiments correlated by a common normalization error, e.g. two detectors observing events from the same beam.

It is straight forward to apply the method sketched in this section also in more complicated situations. For example, if there are several sources of theoretical errors leading to more complicated correlations the pull approach can also be applied by introducing a parameter $`\xi _a`$ for each theoretical error Fogli:2002pt . In a similar way one can treat the case when the compatibility of several experiments should be tested, which are correlated by common theoretical uncertainties. (Consider e.g., the various solar neutrino experiments, which are correlated due to the uncertainties on the solar neutrino flux predictions.)

## VI Conclusions

In this note we have discussed a goodness-of-fit method which was proposed in Ref. Maltoni:2002xd . The so-called parameter goodness-of-fit (PG) can be applied when the global data consists of several statistically independent subsets. Its construction and application are very similar to the standard goodness-of-fit. We gave a formal derivation of the probability distribution function of the proposed test statistic, based on standard theorems of statistics, and illustrated the application of the PG on some examples. We have shown that in the simple case of two data sets determining the mean of a Gaussian, the PG is identical to the intuitive method of testing whether the difference of the two measurements is consistent with zero. Furthermore, we have compared the standard goodness-of-fit and the PG by using real data from neutrino oscillation experiments, which have been the original motivation for the PG. In addition we have illustrated that the so-called pull approach allows to apply the PG also in cases where the data sets are correlated due to theoretical uncertainties.

The proposed method tests the compatibility of different data sets, and it gives sensible results even in cases where the errors are estimated very conservatively and/or the total number of data points is very large. In particular, it avoids the problem that a possible disagreement between data sets becomes diluted by data points which are insensitive to the problem in the fit. The PG can also be very useful when a set consisting of a rather small number of data points is combined with a very large data sample.

To conclude, we believe that physicists should keep an open mind when choosing a statistical method for analyzing experimental data. In many cases much more information can be extracted from data if the optimal statistical tool is used. We think that the method discussed in this note may be useful in several fields of physics, especially where global analyses of large amount of data are performed.

## Acknowledgment

We thank C. Giunti and P. Huber for very useful discussions. Furthermore, we would like to thank many participants of the NOON 2003 workshop in Kanazawa, Japan, for various comments on the PG. M.M. is supported by the Spanish grant BFM2002-00345, by the European Commission RTN network HPRN-CT-2000-00148, by the ESF Neutrino Astrophysics Network N. 86, and by the Marie Curie contract HPMF-CT-2000-01008. T.S. is supported by the “Sonderforschungsbereich 375-95 für Astro-Teilchenphysik” der Deutschen Forschungsgemeinschaft.


text:
# Search for delays between unlike particle emissions in relativistic heavy-ion collisions

## Abstract

The possibilities of unlike particle correlations for a study of particle delays and spatial shifts in relativistic heavy ion collisions are demonstrated. This report represents the unpublished Ref.  in the paper by S. Voloshin et al., Phys. Rev. Lett. 79 (1977) 4766; the report abstract was submitted to the conference Quark Matter’97.

One of the major purposes in studying heavy-ion collisions is to understand the physics of particle interactions at both high particle and high energy densities and its dependence on the colliding energy. Thus, based on the free space cross section (or mean-free-path) arguments, it was suggested long time ago  that the freeze-out sequence is: photons, K+’s, protons and pions. The transport model calculations then confirmed this expectation except for protons which appear to be emitted latest in the mean. This is demonstrated in Fig. 1 where we present the particle yields as functions of the emission time in the reaction c.m.s. for $`Pb+Pb`$ central collisions at 158 A$``$GeV simulated with the RQMD (v2.3) code. Nowadays, in connection with the experiments at SPS with heaviest nuclei and the future experiments at RHIC and LHC, there are of special interest the effects of the eventual phase transitions, like the production of quark-gluon plasma or strangeness distillation, leading to specific delays in particle production. Therefore the measurement of the freeze-out sequence may be able to provide crucial information on the properties of the fireball produced in heavy-ion collisions.

In this report, we address the problem of relative delays and spatial shifts among unlike particles. We use the transport RQMD (v2.3) code to generate the 158 A$``$GeV Pb+Pb central collisions. The correlation functions for unlike particle pairs: ($`\pi ^+,K^+`$), ($`\pi ^+,p`$) and ($`K^+,p`$) are then calculated with the final state interaction code .

In order to minimize the effect of a fast longitudinal motion of the particle sources, we consider here only the mid–rapidity region $`|y|<1`$ and calculate the two–particle correlation functions in the longitudinally co-moving system (LCMS). In this system the particle pairs are emitted transverse ($`x=out`$) to the reaction axis ($`z=long`$). In Fig.2 we present the distributions of the relative $`x`$-coordinates of the emission points and of the relative emission times calculated in RQMD. The distributions of the relative $`y`$\- and $`z`$-coordinates show practically no asymmetry as expected from symmetry reasons (the azimuthal symmetry and the symmetry of the initial system together with the choice of a symmetric mid–rapidity window).

As pointed out in , we can access the information about these asymmetries exploiting the fact that the final state interaction, determining the correlation function of two non–identical particles, depends on the orientation of the relative 4–coordinates of the emission points through the scalar product $`\stackrel{}{q}\stackrel{}{r}^{}`$, where $`\stackrel{}{r}^{}=\stackrel{}{r}_1^{}\stackrel{}{r}_2^{}\{\mathrm{\Delta }x^{},\mathrm{\Delta }y^{},\mathrm{\Delta }z^{}\}`$ and $`\stackrel{}{q}=\stackrel{}{p}_1^{}=\stackrel{}{p}_2^{}`$ are the relative 3–coordinates of the emission points and half the relative momentum in the two–particle c.m.s. Since this system moves in the x-direction with the velocity $`v_{}=|\stackrel{}{p}_1+\stackrel{}{p}_1|/(E_1+E_2)`$ and the Lorentz factor $`\gamma _{}=(1v_{}^2)^{1/2}`$, we have:

$$\mathrm{\Delta }x^{}=\gamma _{}(\mathrm{\Delta }xv_{}\mathrm{\Delta }t),\mathrm{\Delta }y^{}=\mathrm{\Delta }y,\mathrm{\Delta }z^{}=\mathrm{\Delta }z.$$
(1)
Since further $`\stackrel{}{q}\stackrel{}{r}^{}=q_x\mathrm{\Delta }x^{}+q_y\mathrm{\Delta }y^{}+q_z\mathrm{\Delta }z^{}`$ and, in our case $`\mathrm{\Delta }y=\mathrm{\Delta }z=0`$, we can search for the LCMS space (out) and time asymmetries with the help of the correlation functions $`_+`$ and $`_{}`$ corresponding to $`q_x=\stackrel{}{q}\stackrel{}{v}_{}/v_{}>0`$ and $`<0`$, respectively . In the case when the time differences dominate over the spatial ones, i.e. when $`v|t|r`$ on average, the ratio of the correlation functions $`_+`$ and $`_{}`$ directly measures the mean relative emission time, including its sign . Noting that for particles of equal masses the sign of the scalar product $`\stackrel{}{q}\stackrel{}{v}`$ coincides with the sign of the velocity difference $`v_1v_2`$, we can see the simple classical meaning of the above selection. It corresponds to the intuitive expectation of different particle interaction in the case when the faster particle is emitted earlier as compared to the case of its later emission. In former case the interaction between the two particles will be weaker and the correlation function $`_+`$ will be closer to unity than $`_{}`$.

To clarify the origin of the asymmetry effect, we will follow the arguments given in . On the usual assumptions of a small particle phase-space density and sufficiently smooth behavior of the single-particle spectra, the two–particle correlation function is determined by the modulus squared of the two-particle amplitude $`\psi _\stackrel{}{q}(\stackrel{}{r}^{})`$ averaged over the relative coordinates of the emission points in the two-particle c.m.s. :

$$R(p_1,p_2)=|\psi _\stackrel{}{q}(\stackrel{}{r}^{})|^2.$$
(2)
For charged particles the final state interaction is often dominated by the Coulomb interaction so that

$$\psi _\stackrel{}{q}(\stackrel{}{r}^{})=\mathrm{e}^{i\delta _0^c}\sqrt{A_c(\eta )}\mathrm{e}^{i\stackrel{}{q}\stackrel{}{r}^{}}F(i\eta ,1,i\xi ),$$
(3)
where $`\xi =qr^{}+\stackrel{}{q}\stackrel{}{r}^{}qr^{}(1+\mathrm{cos}\theta ^{})`$, $`\eta =1/(k^{}a)`$, $`a=1/(\mu z_1z_2e^2)`$ is the Bohr radius ($`\mu `$ is the reduced mass of the two particles, $`z_1z_2e^2`$ is the product of the particle electric charges; $`a=`$ 248.5, 222.5 and 83.6 fm for the $`\pi ^+p`$, $`K^+p`$ and $`\pi +K^+`$ systems, respectively), $`\delta _0^c=\mathrm{arg}\mathrm{\Gamma }(1+i\eta )`$ is the Coulomb s–wave phase shift,

$$A_c(\eta )=2\pi \eta [\mathrm{exp}(2\pi \eta )1]^1$$
(4)
is the Coulomb factor - modulus squared of the non-relativistic Coulomb wave function at zero distance - introduced already by Fermi in his theory of $`\beta `$–decay,

$$F(\alpha ,1,z)=1+\alpha z+\alpha (\alpha +1)z^2/2!^2+\alpha (\alpha +1)(\alpha +2)z^3/3!^2+\mathrm{}$$
(5)
is the confluent hypergeometric function. At small relative momenta $`qr^{}<1`$ and $`r^{}|a|`$, we can then write

$`R(p_1,p_2)`$ $`=`$ $`A_c(\eta )\left[1+2r^{}(1+\mathrm{cos}\theta ^{})/a+\mathrm{}\right]`$ (6)
$``$ $`A_c(\eta )\left[1+2r^{}/a+2\mathrm{cos}\psi \mathrm{\Delta }x^{}/a+2\mathrm{sin}\psi (\mathrm{sin}\varphi \mathrm{\Delta }y^{}+\mathrm{cos}\varphi \mathrm{\Delta }z^{})/a+\mathrm{}\right],`$ (7)
where $`\psi `$ and $`\varphi `$ are the polar and azimuthal angles of the vector $`\stackrel{}{q}`$ with respect to the velocity vector ($`x`$-axis): $`\stackrel{}{q}=q\{\mathrm{cos}\psi ,\mathrm{sin}\psi \mathrm{sin}\varphi ,\mathrm{sin}\psi \mathrm{cos}\varphi \}`$ and $`\mathrm{cos}\psi =\stackrel{}{q}\stackrel{}{v}_{}/(qv_{})`$. Noting that for particles with the same (opposite) charges the Bohr ”radius” $`a>0`$ ($`<0`$) and the Coulomb factor $`A_c(\eta )`$ approaches unity from below (above), we can immediately see that the interaction (and thus the correlation) is stronger if $`\stackrel{}{q}\stackrel{}{r}^{}=qr^{}\mathrm{cos}\theta ^{}<0`$. In the classical limit of $`qr^{}1`$ this corresponds to the expectation of stronger interaction between particles moving after the emission (on average) towards each other in their c.m.s.*** It should be noted that the two–particle amplitude $`\psi _\stackrel{}{q}(\stackrel{}{r}^{})`$ is often misidentified with the usual two-particle wave function in the scattering problem $`\psi _\stackrel{}{q}(\stackrel{}{r}^{})`$ (see, e.g., ). This misidentification would obviously lead to the opposite asymmetry effect in the correlation function.

Averaging in Eq. (6) over nearly isotropic distribution of the vector $`\stackrel{}{q}`$ at $`q0`$ (nearly uniform $`\mathrm{cos}\psi `$\- and $`\varphi `$–distributions) or recalling that in the considered case $`\mathrm{\Delta }y^{}=\mathrm{\Delta }z^{}=0`$, we get for the ratio $`_+/_{}`$ at $`q0`$

$$_+/_{}1+2\left(\mathrm{cos}\psi _+\mathrm{cos}\psi _{}\right)\mathrm{\Delta }x^{}/a,$$
(8)
where the $`\pm `$ subscripts correspond to the averaging over positive/negative values of $`q_x=q\mathrm{cos}\psi `$. The result of this averaging depends on the selected kinematic region. For example, if only pairs with the relative momentum $`\stackrel{}{q}`$ parallel or antiparallel to the pair velocity are selected, then $`\mathrm{cos}\psi _\pm =\pm 1`$ and the measured asymmetry would be maximal . In the case of averaging over a full range of the uniform $`\mathrm{cos}\psi `$-distribution, we have $`\mathrm{cos}\psi _\pm =\pm 1/2`$ and twice less the asymmetry. For our kinematical selection, the corresponding angular factor in brackets in Eq. (8) is somewhat higher than unity.

Regarding the effect of the strong final state interaction, since the magnitude of the two-particle scattering amplitude $`f`$ is usually much less than the mean $`r^{}`$, it only slightly modifies the correlation functions $`_+`$, $`_{}`$ and their ratio at small $`q`$. In particular, Eq. (8) for the ratio $`_+/_{}`$ at $`q0`$ is modified by the substitution $`\mathrm{\Delta }x^{}\mathrm{\Delta }x^{}+\mathrm{Re}f\mathrm{\Delta }x^{}/r^{}`$. This means that for neutral particles (when $`|a|\mathrm{}`$) the ratio $`_+/_{}1`$ at $`q0`$.In  it was not clearly stated that the asymmetry vanishes at $`q0`$ only on the absence of the Coulomb interaction.

The correlation functions $`_+`$, $`_{}`$ and their ratio are plotted in Fig. 3. We can see that for $`\pi ^+p`$ and $`\pi ^+K^+`$ systems these ratios are less than unity at small values of $`q`$, while for $`K^+p`$ system the ratio $`_+/_{}`$ practically coincides with unity. These results well agree with the mean values of $`\mathrm{\Delta }t`$, $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }x^{}`$ presented in Table 1. We may see from Eqs. (1) and (8) that the absence of the effect in the $`_+/_{}`$ ratio for the $`K^+p`$ system is due to practically the complete compensation of the space and time asymmetries leading to $`\mathrm{\Delta }x^{}0`$. For $`\pi ^+p`$ system the effect is determined mainly by the x-asymmetry. For $`\pi ^+K^+`$ system both the x- and time-asymmetries contribute in the same direction, the latter contribution being somewhat larger.

We have thus demonstrated that the unlike particle correlations allow to directly access the information about space-time asymmetries of particle production in ultrarelativistic heavy ion collisions even at present energies. The separation of the relative time delays from the spatial asymmetry is, in principle, possible (see Eq. (1)) by studying the ratio $`_+/_{}`$ in different intervals of the pair velocity. However, such a study can hardly be done in a model independent way due to the velocity dependence of the out asymmetry $`\mathrm{\Delta }x`$ (particularly, both $`\mathrm{\Delta }x`$ and $`v_{}\mathrm{\Delta }t`$ vanish at zero LCMS pair velocity $`v_{}`$).