text:
# Solar Neutrinos Before and After KamLAND

## 1 Introduction

KamLAND is the first experiment to explore with a terrestrial beam the region of neutrino oscillation parameters that is relevant for solar neutrino oscillations. This spectacular achievement  makes possible new studies of the physics of neutrinos  and of the solar interior . We concentrate here on refining the allowed regions in neutrino oscillation space and on what can be learned about the total 8B solar neutrino flux, as well as the sterile neutrino component of the 8B neutrino flux. The procedures we use have been described in refs. .

We assume throughout this paper that CPT is satisfied and that the anti-neutrino measurements by KamLAND apply directly to the neutrino section. In fact, the first KamLAND results already show that CPT is satisfied in the weak sector to a characteristic accuracy of about $`10^{20}`$ GeV .

We begin in section 2 by determining the allowed oscillation regions using the recently released KamLAND data. Our principal results, which are in good agreement with the analysis of the KamLAND collaboration , are summarized in figure 1 and in table 3. We also show that the numerical values of the physical quantities that we calculate in this paper are insensitive to the details of the analysis of the KamLAND data (see discussion at the end of section 6).

We derive in section 3 a global solution for the neutrino oscillation parameters using, together with the KamLAND measurements, all the available solar neutrino data from the chlorine , gallium , Super-Kamiokande , and SNO  experiments. The enormous reduction in the size of the allowed oscillation regions is evident from a comparison of the ‘Before’ and ‘After’ allowed regions, figure 3 and figure 4.

We determine in section 4 the total 8B solar neutrino flux, as well as the uncertainties in the flux that are implied by the existing experimental data. This flux can be compared directly with the flux predicted by the standard solar model.

Could there be a large but unobservable flux of 8B solar neutrinos that reaches Earth in the form of sterile (right-handed) neutrinos? Prior to the announcement of the KamLAND results , the answer to this question was a resounding “Yes”. At $`3\sigma `$, the sterile component of the 8B flux could be, prior to KamLAND, as large as 1.1 times the BP00 predicted flux, which amounts to 50% of the total flux for this extreme case. However, the KamLAND measurement allows us to place a strong upper limit on the flux of sterile 8B neutrinos. We determine in section 5 the best-estimate sterile 8B neutrino flux and the associated uncertainties in the sterile flux.

We summarize and discuss our main conclusions in section 6.

## 2 KamLAND only

Figure 1 shows the allowed regions for anti-neutrino oscillations that we have determined using the description of the detector operation and the measured data presented in the recent KamLAND publication . The left hand panel of figure 1 shows the allowed regions calculated using only the KamLAND total event rate; the right hand panel shows the allowed regions calculated using the experimental energy spectrum, which also contains a rate normalization. In computing the right hand panel of figure 1, we also used the data from the CHOOZ reactor experiment. This experiment rules out solutions with $`\mathrm{\Delta }m^2`$ greater than $`1\times 10^3\mathrm{eV}^2`$ at $`3\sigma `$ or $`8\times 10^4\mathrm{eV}^2`$ at 99% CL.

Section 2.1 contains our version (used in the following sections of the paper) of the science results from the KamLAND measurements and may be of general interest. Section 2.2 contains a description of how we do the data analysis and is rather specialized. We describe in sec. 2.3 Earth matter effects in the KamLAND experiment. The results of sec. 2.3 show that matter effects must be included in future precision analyses of the KamLAND experiment. We have carried out the KamLAND-only and the solar plus KamLAND analyses described in this paper both with and without including matter effects for the KamLAND experiment and have verified that, at the present level of accuracy for KamLAND measurements, the predictions for future experimental measurements are unaffected by including matter effects for KamLAND. Section 2.2 and Section 2.3 are rather technical and can be skipped by the reader who is primarily interested in the science.

We have verified that different plausible approaches to analyzing the KamLAND data do not lead to significant differences in the quantities in which we are most interested: the total 8B solar neutrino flux, the sterile component of the 8B neutrino flux, the CL at which maximal mixing is excluded, and the predictions for future solar neutrino experiments.

### 2.1 The KamLAND allowed regions

As expected, the allowed regions shown in figure 1 are in agreement with those reported by the KamLAND collaboration. It was necessary to establish explicitly this agreement in order to have confidence in the subsequent steps of our analysis, which involve determining: the globally allowed regions (section 3.1.1), predictions for future solar neutrino experiments (section 3.1.2), the total 8B solar neutrino flux (section 4), an upper limit to the sterile component of the 8B neutrino flux (section 5), and an upper limit to the CNO contribution of the solar luminosity (ref. ).

The results we obtain for the global solar plus KamLAND contours of the allowed regions, and all of the quantities derived from these contours, are essentially independent of the details of the analysis of the KamLAND results. We have verified this independence by analyzing the KamLAND measurements in a number of different ways (see discussion at the end of section 6).

The most important aspect of figure 1 is the demonstration by the KamLAND results that anti-neutrinos oscillate with parameters that are consistent with the LMA solar neutrino solution . We shall exploit this extraordinary achievement in detail in the following sections. The symmetry shown in figure 1 with respect to the line $`\theta =\pi /4`$ (or $`\mathrm{tan}^2\theta =1`$) is just a reflection of the fact that the KamLAND experiment is an essentially vacuum experiment (but cf. sec. 2.3 and fig. 2) and hence the survival probabilities depend almost entirely upon $`\mathrm{sin}^22\theta `$.

For the rates-only analysis, the best-fit oscillation parameters form a degenerate line in the $`\mathrm{\Delta }m^2`$-$`\mathrm{tan}^2\theta `$ plane. For large $`\mathrm{\Delta }m^2`$, the line only depends upon $`\theta `$ and is described by $`\mathrm{tan}^2\theta =0.36`$.

Expressed in terms of the rate expected if there were no flavor changes, the KamLAND rate is

$$R_{\mathrm{KamLAND}}=0.611\pm 0.094.$$
(1)
The measured rate is in excellent agreement with the pre-KamLAND expectations

$$R_{\mathrm{KamLAND},\mathrm{expected}}=0.58_{0.27}^{+0.10},$$
(2)
based on purely solar neutrino experiments.

As can be seen in the right hand panel of figure 1 and as is summarized in table 1, the KamLAND spectral data lead to three local minima. The allowed region is separated into ‘islands.’ These islands correspond to oscillations with wavelengths that are approximately tuned to the average distance between the reactors and the detector, $`180`$ km. The best-fit point in the lowest mass island ($`\mathrm{\Delta }m^2=1.5\times 10^5\mathrm{eV}^2`$) is near the first-maximum in the oscillation probability (minimum in the event rate). The overall best-fit point ($`\mathrm{\Delta }m^2=7.1\times 10^5\mathrm{eV}^2`$) lies within the island around the second maximum in the oscillation probability. The tuning is less accurate for the best-fit point in the highest mass island.

For each of the allowed islands in neutrino oscillation parameter space that are shown in the right hand panel of figure 1 (analysis of the energy spectrum), table 1 gives the best-fit values of $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$. The global best-fit point for the KamLAND data set lies at $`\mathrm{\Delta }m^2=7.1\times 10^5`$ eV2 and mixing angle $`\mathrm{tan}^2\theta =0.52`$. For the same mass difference, maximal mixing is slightly less favored, $`\mathrm{\Delta }\chi ^2=0.4`$. With the present statistical accuracy, KamLAND cannot discriminate between large angle (non-maximal) mixing and maximal mixing.

At what level are other, non-LMA oscillation solutions excluded? For $`\mathrm{\Delta }m^2<10^6\mathrm{eV}^2`$, there is no significant oscillation probability in the KamLAND experiment. All such low $`\mathrm{\Delta }m^2`$ solutions (like the previously-discussed LOW and Vacuum solar neutrino oscillation solutions), have $`\mathrm{\Delta }\chi ^2>15.9`$ and are excluded by at least $`3.6\sigma `$ (2 dof).

Comparing the total number of observed and expected events in KamLAND we find that alternative explanations to the solar neutrino problem that predict no deficit in the KamLAND experiment are now excluded at 3.6$`\sigma `$. This applies, for instance, to spin flavor precession, or flavor changing neutrino interactions which are now ruled out at 3.6$`\sigma `$ as the dominant mechanism for solar neutrino flavor conversion. In particular, all ‘non-standard solar model’ explanations of the ‘solar neutrino problem’ which assume particle physics with $`P_{ee}=1`$ (no neutrino oscillations) are also excluded by the KamLAND experiment at $`3.6\sigma `$.

### 2.2 KamLAND analysis

¿From the experimental point of view, KamLAND involves different systematic uncertainties than solar neutrino experiments. From a theoretical point of view, KamLAND uses anti-neutrinos to perform a pure disappearance experiment in vacuum. Solar neutrino studies use neutrinos to perform disappearance (or appearance) experiments in matter (or in vacuum). The analysis of the KamLAND data requires a treatment that is different from, and in many ways more simple than, the global analysis of solar neutrino data.

Using the data provided in the KamLAND paper , we calculate the positron spectrum in KamLAND detector with the procedures described in refs. . In the absence of neutrino oscillations, we find (in agreement with ref. ) 86.8 expected neutrino events above 2.6 MeV visible energy, for 145.1 days of live time, a 408 ton fiducial mass with 3.46 $`\times 10^{31}`$ free target protons, detection efficiency of 78.3%, energy resolution $`\sigma (E)/E=7.5\%/\sqrt{E}`$, integrated total thermal power flux during the measurement live time of 254 Joule/$`\mathrm{cm}^2`$, and relative fission yields from different fuel components as specified in ref. . The positron energy spectrum that we calculate is in excellent agreement with the energy spectrum without oscillations presented by the KamLAND collaboration.

We perform, following the KamLAND collaboration, two analyses of the first KamLAND results: 1) including only the total event rate, and 2) including the total measured energy spectrum. For the rate only analysis, the $`\chi ^2`$ can be written as ,

$$\chi _{\mathrm{Rate}}^2=\frac{(0.611R(\theta ,\mathrm{\Delta }m^2))^2}{\sigma _\mathrm{R}^2},$$
(3)
with $`\sigma _\mathrm{R}^2=(0.085)^2+(0.041)^2`$ and

$$R=1\mathrm{sin}^22\theta \mathrm{sin}^2\left(\frac{\mathrm{\Delta }m^2L}{4E_\nu }\right).$$
(4)
We denote by $``$ the convolution over the energy spectrum, the interaction cross section, and the energy resolution function. For simplicity in the notation, we have neglected in writing eq. 4 all matter effects (see sec. 2.3 for details of the matter effects).

In the left hand panel of figure 1, we show the allowed regions (2 dof) that are obtained from this rate-only analysis. The vertical asymptotes can be calculated simply from eq. (3) since for large $`\mathrm{\Delta }m^2`$ the quantity $`\mathrm{sin}^2\left(\mathrm{\Delta }m^2L/4E_\nu \right)`$ averages to 0.5 . Thus in the large $`\mathrm{\Delta }m^2`$ limit, the mixing angle boundary of the allowed region satisfies $`\mathrm{sin}^2(2\theta )=2\sigma _\mathrm{R}\sqrt{\mathrm{\Delta }\chi _{CL}^2}+0.778`$. Here $`\mathrm{\Delta }\chi ^2`$ is fixed by the CL of interest; $`\mathrm{\Delta }\chi ^2=5.99(11.83)`$ for 95% CL (for $`3\sigma `$ CL).

We included the data from the KamLAND energy spectrum in a way that has become conventional when analyzing low statistics binned data , assuming that data is Poisson distributed. We define the $`\chi ^2`$ function as

$`\chi ^2(\mathrm{\Delta }m^2,\theta )=2=`$ (5)
$`\mathrm{min}_\alpha {\displaystyle \underset{i}{}}\left[2(\alpha N_i^{\mathrm{th}}(\mathrm{\Delta }m^2,\theta )N_i^{\mathrm{exp}})+2N_i^{\mathrm{exp}}\mathrm{ln}\left({\displaystyle \frac{N_i^{\mathrm{exp}}}{\alpha N_i^{\mathrm{th}}(\mathrm{\Delta }m^2,\theta )}}\right)\right]+{\displaystyle \frac{(\alpha 1)^2}{\sigma _{\mathrm{sys}}^2}},`$
where $`\alpha `$ is an absolute normalization constant and $`\sigma _{\mathrm{sys}}=6.42\%`$ represents the systematic uncertainties presented in table 2 of ref. . The KamLAND collaboration analyzed the event sample by separating the total rate from the spectrum shape and using a maximum likelihood method.

In the right hand panel of figure 1, we show the allowed regions (2 dof) that we obtain from our analysis of the KamLAND spectrum.

We have performed a series of tests of the dependence of the allowed solutions on the statistical treatment of the data. For example, we have performed a standard $`\chi ^2`$ analysis using the first 5 energy bins above 2.6 MeV visible energy as presented by the KamLAND collaboration, but with all the events above the fifth bin ($`E>4.725`$ MeV) collected into a single bin (or into two bins). We then assumed Gaussian statistics for all bins. We find that the exact size and shape of the $`3\sigma `$ allowed regions of the KamLAND-only analysis depend mildly on the strategy adopted, i.e., whether we use eq. (5) or the standard $`\chi ^2`$. However, the differences are unimportant for quantities calculated using the global solar plus KamLAND allowed regions. We have also studied the effect of carefully treating the energy dependence of the energy threshold error by computing, as a function of the oscillation parameters, the systematic uncertainty for each bin due to the estimated systematic error for the energy scale (1.91% at 2.6 MeV). This refinement has no discernible effect on any of the quantities we calculate in this paper.

We discuss at the end of section 6 the small quantitative dependence upon analysis strategy of each of the quantities we calculate in this paper.

### 2.3 Matter effects in the KamLAND experiment

It is conventional to ignore matter effects in analyzing the KamLAND results. However, as we shall see below, matter effects can influence the calculated anti-neutrino event rate by as much as $`4`$% if the mass and mixing angle are near the region of maximum sensitivity. Thus matter effects will have to be included in future analyses of the KamLAND results when greater statistical and systematic precision are obtained.

In the results described in this paper, we have included matter effects by assuming that the anti-neutrinos traverse a constant density medium ($`\rho =2.7\mathrm{gm}\mathrm{cm}^3`$, cf. ref. )with the path length appropriate to each reactor. Since matter effects are small, the assumption of a constant density is an excellent approximation. Once matter effects are included, the survival probability for the KamLAND experiment also depends on the active-sterile admixture. However, allowing sterile neutrinos reduces the influence of matter effects with respect to what is calculated for the pure active case.

To the accuracy that can be discerned by even the most expert eye, matter effects do not change the contours of figures given in the present paper.

We define, $`F(\mathrm{matter}\mathrm{vs}\mathrm{vacuum})`$, as the fractional difference in the event rate for the KamLAND detector calculated with and without including matter effects in the Earth. Thus

$$F(\mathrm{matter}\mathrm{vs}\mathrm{vacuum})=\frac{R_{\mathrm{matter}}R_{\mathrm{vacuum}}}{R_{\mathrm{vacuum}}}.$$
(6)
Figure 2 shows the contours of $`F(\mathrm{matter}\mathrm{vs}\mathrm{vacuum})`$ in the $`\mathrm{\Delta }m^2`$-$`\mathrm{tan}^2\theta `$ plane for active neutrinos. The maximum value of $`|F(\mathrm{matter}\mathrm{vs}\mathrm{vacuum})|`$ corresponds to -4.2% (+4.3%) at $`\mathrm{\Delta }m^2=3.55\times 10^5\mathrm{eV}^2`$ and $`\mathrm{tan}^2\theta =0.43(2.18)`$. Within the $`1\sigma (3\sigma )`$ allowed region of the analysis of the KamLAND energy spectrum (cf. the right hand panel of fig. 1), the maximum value of $`|F(\mathrm{matter}\mathrm{vs}\mathrm{vacuum})|`$ corresponds to $`0.9`$% ($`3.9`$%). The maximum change in $`\chi ^2`$ due to including matter effects in the KamLAND-only analysis is 0.2 (1.5) at $`1\sigma (3\sigma )`$. At the best-fit point for the analysis of the KamLAND energy spectrum, matter effects change the value of $`\chi ^2`$ by only 0.1%.

We anticipate the results given latter in sec. 3 (see especially fig. 4) by noting that the matter effects in the KamLAND experiment are somewhat less important within the globally allowed solar plus KamLAND oscillation regions. For the global solution, $`|F(\mathrm{matter}\mathrm{vs}\mathrm{vacuum})|`$ varies between $`0.9`$% ($`2.7`$%) at $`1\sigma (3\sigma )`$. The maximum change in $`\chi ^2`$ in the global analysis due to matter effects is 0.2 (1.0).

For simplicity, we have described in the rest of this paper the KamLAND experiment as if matter effects could be completely ignored.

## 3 Global solutions: solar plus KamLAND

We present and discuss in section 3.1 the allowed regions in neutrino oscillation parameter space that are permitted at $`3\sigma `$ by the combined solar and KamLAND experimental data (see especially figure 4). We also present in this section (see especially table 2) the predictions of future solar neutrino observables that are calculated using the current global oscillation solutions.

In order to isolate the crucial measurements that have identified LMA as the correct description of solar neutrino oscillations, we present in section 3.2 the results of an analysis of just the rates of the solar neutrino experiments. We show in this section that, if we did not have the essentially undistorted electron recoil energy spectrum measured by Super-Kamiokande and confirmed by SNO, vacuum oscillations would be the preferred solution prior to the KamLAND measurements.

We include two sections that are primarily intended for aficionados of neutrino oscillation analyses. We outline briefly in section 3.3 how we have combined the solar and the KamLAND analyses. We have already described in section 2.2 how we have analyzed the KamLAND measurements. In section 3.3.1, we summarize our treatment of the solar data. We describe in section 3.3.2 our treatment of the sterile component of the solar neutrino flux.

### 3.1 Scientific results

We describe in section 3.1.1 the allowed regions in neutrino oscillation space that are determined by the solar and KamLAND data and in section 3.1.2 we present predictions for future solar neutrino experiments of the presently allowed global oscillation solutions.

#### 3.1.1 Allowed regions

Figure 3 presents the ‘Before’ figure, the allowed neutrino oscillation regions for all the available solar neutrino data prior to the announcement of the first KamLAND results. We show in figure 3 the results of our analysis for this paper of the pre-KamLAND solar neutrino data. We include here for the first time the recently reported GNO measurements ; this inclusion does not affect figure 3 significantly. Our calculational procedures, which are described in ref. , contain some refinements not taken account of in other analyses. However, the overall results for the pre-KamLAND allowed regions are very similar in most treatments, cf. .

Figure 4 displays the ‘After’ figure, the much reduced allowed regions that are determined by including the recently announced KamLAND measurements together with the solar neutrino data.

Comparing figure 3 with figure 4, we see that KamLAND has shrunk enormously the allowed parameter space that is accessible to theoretical models, by five orders of magnitude in $`\mathrm{\Delta }m^2`$ and one order of magnitude in $`\mathrm{tan}^2\theta `$. The LMA region in parameter space is the only remaining allowed solution. The best-fit LOW solution is excluded at $`4.8\sigma `$ and vacuum solutions are excluded at $`4.9\sigma `$. The once ‘favored’ SMA solution is now excluded at $`6.1\sigma `$.

The allowed range for the mixing angle within the 1$`\sigma `$ (3$`\sigma `$) region shown in the ‘After’ figure 4 is

$$0.36(\mathrm{\hspace{0.17em}0.27})<\mathrm{tan}^2\theta <0.58(0.94).$$
(7)
For the pre-KamLAND analysis, figure 3, the allowed $`3\sigma `$ domain for the mixing angle was $`0.26<\mathrm{tan}^2\theta <0.98`$. Therefore, the allowed range of mixing angles has only been marginally improved by the inclusion of the KamLAND data.

At what confidence level is maximal mixing excluded? There are no solutions with $`\mathrm{tan}^2\theta =1`$ at the $`3.1\sigma `$ CL. For the pre-KamLAND analysis, figure 3, there were no solutions with maximal mixing within the LMA region at $`3.1\sigma `$. However, maximal mixing was still allowed at $`2.8\sigma `$ within the vacuum solution regions, which is now excluded.

Figure 5 shows how solar neutrino experiments break the symmetry of the survival probability, $`P_{\mathrm{ee}}\left(E\right)`$, about a mixing angle of $`\theta =\pi /2`$ that exists for vacuum oscillation experiments. Vacuum oscillation survival probabilities are identical for two mixing angles $`\theta _1`$ and $`\theta _2`$, with $`\theta _2=\pi /2\theta _1`$ and a fixed value for $`\mathrm{\Delta }m^2`$. The left hand panel shows the vacuum survival probabilities for a source-detector separation of 180 km, similar to the actual experimental situation for KamLAND. The right hand panels displays the daytime survival probabilities for the two LMA solar neutrino oscillation solutions. Matter effects in the Sun are responsible for the different LMA survival probabilities.

Which solar neutrino experiments choose the lower probabilities (the smaller value for $`\mathrm{tan}^2\theta `$)? Originally, the chlorine experiment, in conjunction with the standard solar model predictions, indicated strongly that the smaller probabilities were correct. The measured event rate in the chlorine experiment, $`2.6`$ SNU , was closer to $`0.35`$ of the standard solar model prediction than to $`0.65`$ of the solar model prediction ($`8`$ SNU) . This result was confirmed by the comparison of the Super-Kamiokande $`\nu e`$ scattering measurements with the SNO CC measurement . The comparison within the same experiment of the SNO CC and NC rates selects unambiguously the curve with the lower survival probabilities .

Finally, we have evaluated the present status of ‘non-standard’ solar model explanations of the solar plus KamLAND data under the hypothesis of standard particle physics. We perform a fit to the total measured rates in the chlorine , gallium , Super-Kamiokande , SNO CC , and SNO NC  experiments, allowing the dominant $`pp`$, 8B, 7B and CNO fluxes to vary freely but imposing the luminosity constraint (see  and references therein for details). We find that $`\chi _{\mathrm{min},\mathrm{solar}}^2=53.9`$ for $`2`$ dof (5 rates - 4 free fluxes + 1 constraint), which implies that non-standard solar model explanations of the solar data can be excluded at $`7\sigma `$. Adding to the analysis the evidence of antineutrino disappearance observed in the KamLAND experiment, we find that non-standard solar models invented to avoid solar neutrino oscillations are excluded by solar plus KamLAND experiments at $`7.9\sigma `$.

#### 3.1.2 Predictions for future experiments

Table 2 presents the best-fit predictions and the expected ranges for some important solar neutrino observables. The results summarized in table 2 correspond to the globally allowed oscillation regions shown in figure 4. We consider predictions for SNO, for a large SK-like neutrino-electron scattering experiment, and for a neutrino-electron scattering detector, BOREXINO  or KamLAND , that observes the 7Be solar neutrinos. We also include predictions for a generic $`pp`$ neutrino detector. The notation used here is the same as in ref. .

The day-night asymmetry $`A_{\mathrm{N}\mathrm{D}}`$ is defined in terms of the Day and the Night event rates by the expression

$$A_{\mathrm{N}\mathrm{D}}=2\frac{[\mathrm{Night}\mathrm{Day}]}{[\mathrm{Night}+\mathrm{Day}]}.$$
(8)
The prediction for the day-night asymmetry in the SNO CC measurement is one of the most interesting new results to come from the much smaller, post-KamLAND allowed regions (see discussion of earth matter effects in the LMA regime in refs. ).

Within the unique 1$`\sigma `$ region (2 dof), the prediction is

$$A_{\mathrm{N}\mathrm{D}}(\mathrm{SNO}\mathrm{CC})=3.3_{0.6}^{+0.7}\%.$$
(9)
Within the $`3\sigma `$ region (2 dof), $`A_{\mathrm{N}\mathrm{D}}(\mathrm{SNO}\mathrm{CC})=3.3_{3.1}^{+5.7}\%.`$

An accurate day-night asymmetry measurement in the SNO CC mode could, in principle, help to discriminate between larger and smaller values of $`\mathrm{\Delta }m^2`$. For values of $`\mathrm{\Delta }m^2<10^4\mathrm{eV}^2`$ within the currently allowed $`3\sigma `$ (2 dof) region,

$$A_{\mathrm{N}\mathrm{D}}(\mathrm{SNO}\mathrm{CC})=3.3_{1.6}^{+5.7}\%,(\mathrm{\Delta }m^2<10^4eV^2).$$
(10)
For values of $`\mathrm{\Delta }m^210^4\mathrm{eV}^2`$, the day-night asymmetry is small,

$$A_{\mathrm{N}\mathrm{D}}(\mathrm{SNO}\mathrm{CC})=0.6_{0.4}^{+1.1}\%,(\mathrm{\Delta }m^210^4eV^2).$$
(11)
Unfortunately, it will be very difficult to reach the required experimental precision.

In table 2 we also present the predicted day-night asymmetry at a SK-like neutrino-electron scattering detector with a threshold $`E_{\mathrm{th}}=5`$ MeV. For ES the day-night effect is decreased relative to pure CC interactions by the contribution of the neutral currents.

The prediction for the reduced 7Be $`\nu e`$ scattering rate,

$$[{}_{}{}^{7}\mathrm{Be}]\frac{\mathrm{Observed}\nu \mathrm{e}\mathrm{scattering}\mathrm{rate}}{\mathrm{BP00}\mathrm{predicted}\mathrm{rate}},$$
(12)
is remarkably precise and remarkably stable. Within the $`1\sigma `$ allowed region, the pre-KamLAND prediction was $`[{}_{}{}^{7}\mathrm{Be}]=0.64\pm 0.03,1\sigma `$ . The post-KamLAND prediction for the currently allowed $`1\sigma `$ region is

$$\left[{}_{}{}^{7}\mathrm{Be}\right]_{\mathrm{post}\mathrm{KamLAND}}=0.64\pm 0.02.$$
(13)
The predicted event rate for the 7Be rate experiment has a precision of $`\pm 3`$% in units of the rate calculated with the BP00 flux and assuming no neutrino oscillations.

Now that LMA is the only oscillation solution allowed at $`3\sigma `$, the predicted 7Be day-night asymmetry averaged over a year is too small to be observed. Within the current 3$`\sigma `$ allowed oscillation region, the limit on the predicted values of $`A_{\mathrm{N}\mathrm{D}}`$ is

$$\left|A_{\mathrm{N}\mathrm{D}}({}_{}{}^{7}\mathrm{Be})\right|<0.001.$$
(14)
We also include in table 2 predictions for the rate at which 8B neutrinos will be scattered by electrons resulting in recoil electrons in BOREXINO or KamLAND above two different electron recoil kinetic-energy energy thresholds, $`T_{\mathrm{th}}=3.5`$ MeV and $`T_{\mathrm{th}}=5`$ MeV. The predictions are given in terms of the reduced event rate for 8B neutrinos, \[R(8B)\], which is defined analogously to \[7Be\] \[see eq. (12)\] relative to the scattering rate predicted for the standard solar model flux with no oscillations. The predicted ranges of \[R(8B)\] are obtained from the allowed 1$`\sigma `$ or 3$`\sigma `$ regions in the three parameter space of $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta `$, and $`f_{\mathrm{B},\mathrm{total}}`$. Not surprisingly, the predicted 8B neutrino-electron scattering rates are essentially identical to the neutrino-electron scattering rates measured by the SNO  and Super-Kamiokande  experiments, which are $`0.47\pm 0.05`$ and $`0.465\pm 0.015`$, respectively. For plausible assumptions about how the BOREXINO detector will operate, one expects  $`175`$ events above a $`3.5`$ MeV recoil electron energy threshold and $`110`$ events per year above a $`5`$ MeV threshold.

In addition, we have evaluated the predictions for a generic $`pp`$ neutrino-electron scattering detector. The reduced rate for this detector is defined by the relation

$$[pp]\frac{\mathrm{Observed}\nu \mathrm{e}\mathrm{scattering}\mathrm{rate}}{\mathrm{BP00}\mathrm{predicted}\mathrm{rate}}.$$
(15)
We present the predicted rate for two plausible kinetic energy thresholds, $`100`$ keV and $`50`$ keV. The predicted rate is precise and robust. The precision of the prediction is $`\pm 3`$%; the best-fit predictions have changed by less than $`1`$% from the pre-KamLAND predictions .

The principal differences between the post-KamLAND predictions summarized in table 2 and the pre-KamLAND predictions given in table 2 of ref.  are a consequence of the disappearance of the LOW and the vacuum solutions. As a consequence, no significant day-night asymmetry nor seasonal variations are expected for 7Be neutrino measurements by BOREXINO or KamLAND. For all other observables, the main effect is the reduction of the 1$`\sigma `$ LMA ranges. Within the 3$`\sigma `$ LMA ranges, the changes are generally insignificant.

### 3.2 Rates only analysis

It is instructive to carry out a global analysis using only the total measured rates in the chlorine , gallium , Super-Kamiokande , SNO CC , and SNO NC  experiments. In this simplified analysis, we ignore all spectral energy data and all measurements of the zenith-angle (or day-night) dependence. We let the total 8B neutrino flux be a free variable in the standard rates-only analysis. For simplicity, we consider pure active, or pure sterile, neutrino oscillations. Therefore, the rates-only analysis has three degrees of freedom ($`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta `$, and the 8B neutrino flux) and five measured rates.

The best-fit solution lies in the vacuum region. The oscillation parameters of the best fit solution are: $`\mathrm{\Delta }m^2=8\times 10^{11}\mathrm{eV}^2,\mathrm{tan}^2\theta =0.30(3.35)`$ with $`\chi _{\mathrm{min}}^2=2.2`$. The values of $`\chi _{\mathrm{min}}^2`$ for the other well known solutions are: SMA (3.1) , LMA (3.8), and LOW (8.6), and pure sterile (23.1). All four solutions for active neutrinos, VAC, SMA, LMA, and LOW are allowed at $`3\sigma `$ in the rates-only analysis. If only rates are considered, the VAC, SMA, and LMA solutions are essentially equivalent.

We have done several tests to determine the robustness of this conclusion. The vacuum solutions always appear to be slightly preferred in the rates-only analyses, but the relative order of SMA-LMA can be reversed by slightly changes (such as using the original SNO CC data  with a 6.75 MeV threshold rather than the more recent data  with a 5.0 MeV threshold).

This analysis illustrates the fact that the measurement of the essentially undistorted recoil electron energy spectrum by Super-Kamiokande , which was confirmed by SNO , was required in order to establish, prior to the KamLAND measurements, that LMA was the favored solution. Previous rates-only analyses have established the critical importance of the Super-Kamiokande spectrum analysis . However, this is the first calculation that we know about in which all five of the currently available rate experiments have been included.

### 3.3 KamLAND plus Solar Analysis

It is convenient and conventional to discuss the experimentally determined 8B neutrino flux in units of the flux predicted by the standard solar model. Let $`f_{\mathrm{B},\mathrm{total}}`$ be the ratio of the total 8B solar neutrino flux to the flux predicted by the BP00 standard solar model, i.e.,

$$f_{\mathrm{B},\mathrm{total}}\frac{\varphi (^8\mathrm{B})}{\varphi (^8\mathrm{B})_{\mathrm{BP00}}},$$
(16)
where $`\varphi (^8\mathrm{B})_{\mathrm{BP00}}5.05\times 10^6\mathrm{cm}^2\mathrm{s}^1`$ .

We allow for the possibility that the 8B neutrino flux that arrives at Earth contains a significant component of sterile neutrinos which we parameterize in terms of the active-sterile admixture parameter $`\mathrm{sin}^2\eta `$ (See Sec. 3.3.2).

We calculate the global $`\chi ^2`$ by fitting to all the available data, solar plus reactor. Formally, the global $`\chi ^2`$ can be written in the form

$$\chi _{\mathrm{global}}^2=\chi _{\mathrm{solar}}^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,f_{\mathrm{B},\mathrm{total}},\mathrm{sin}^2\eta )+\chi _{\mathrm{KamLAND}}^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ).$$
(17)
There are four free parameters in $`\chi _{\mathrm{global}}^2`$, $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,f_{\mathrm{B},\mathrm{total}},`$ and $`sin^2\eta `$, although $`\chi _{\mathrm{KamLAND}}^2`$ depends only on $`\mathrm{\Delta }m^2,`$ and $`\mathrm{tan}^2\theta `$. The most computationally intensive task in our analysis of solar plus KamLAND data, computing the multi-dimensional $`\chi ^2`$ for all of the solar neutrino data, was calculated before the KamLAND announcement  and stored. This data set amounted to about 9.4 gigabytes. We do not include the CHOOZ data  explicitly in eq. 17 since the combined solar data also excluded larger values of $`\mathrm{\Delta }m^2`$.

We have carried out a global analysis  of the combined solar neutrino data (80 measurements), see the discussion in section 3.3.1 below, and the KamLAND data (as described in section 2.2). We also consider the possibility that the 8B neutrino flux that arrives at Earth contains a significant component of sterile neutrinos. We describe in section 3.3.2 the aspects of the analysis that are related to sterile neutrinos.

#### 3.3.1 Solar analysis

The details of our calculational methods are described in ref. . We treat the 8B solar neutrino flux as a free parameter, but adopt the predictions of the BP00 solar model  for the other neutrino fluxes and uncertainties. The free 8B strategy adopted in the present paper was, in ref. , the preferred analysis procedure among the three analysis strategies that were used earlier.

We summarize briefly in this subsection the main features of our analysis of solar neutrino data.

For the solar neutrino analysis, we use 80 data points. We include the two measured radiochemical rates, from the chlorine  and the gallium  experiments, the 44 zenith-spectral energy bins of the electron neutrino scattering signal measured by the Super-Kamiokande collaboration , and the 34 day-night spectral energy bins measured with the SNO  detector. The average rate for the gallium experiments  has been updated to $`70.8\pm 4.4`$ SNU by including the preliminary GNO data reported at Neutrino 2002 .

We take account of the BP00  predicted fluxes and uncertainties for all solar neutrino sources except for 8B neutrinos. The Super-Kamiokande  and SNO  experiments can now determine the 8B solar neutrino flux at a level of precision that is comparable with the precision of the BP00 prediction. Therefore, we treat the total 8B solar neutrino flux as a free parameter to be determined by experiment and to be compared with solar model predictions.

Following the SNO collaboration , we have considered together (consistent with the way they are measured) the SNO charged current, neutral current, and electron scattering events. For each point in neutrino oscillation parameter space, we have computed the 34 day-night spectrum energy bins by convolving the computed survival probabilities with the appropriate neutrino cross sections, as well as the experimental energy resolution and sensitivity functions.

We include the errors and their correlations for all of these observables. We construct an $`80\times 80`$ covariance matrix that includes the effect of correlations between the differ error as off diagonal elements. The errors and their correlations are described in in the Appendix of ref.  and in section 5.4 of ref.

The analysis of the theoretical and experimental uncertainties is subtle. Many of the differences in the precise borders of the allowed oscillation regions that are reported by different groups may be traced to variations in the ways in which the errors are treated. We have performed the calculations accurately, although some of the refinements we have used are computationally demanding. As described in the Appendix of ref. , we include the energy dependencies and the correlations of the chlorine and gallium neutrino absorption cross sections. For both the Super-Kamiokande and the SNO experiments, we take account of the energy dependence and the correlation of the 8B spectral energy shape errors, the absolute energy scale errors, and the energy resolution errors.

#### 3.3.2 Sterile component

We have also considered the possibility that $`\nu _e`$ oscillates into a state that is a linear combination of active ($`\nu _a`$) and sterile ($`\nu _s`$) neutrino states,

$$\nu _e\mathrm{cos}\eta \nu _a+\mathrm{sin}\eta \nu _s,$$
(18)
where $`\eta `$ is the parameter that describes the active-sterile admixture. This admixture arises in the framework of 4-$`\nu `$ mixing . The total 8B neutrino flux can be written

$$\varphi ({}_{}{}^{8}\mathrm{B})_{\mathrm{total}}=\varphi (\nu _e)+\varphi (\nu _a)+\varphi (\nu _s),$$
(19)
where $`\varphi (\nu _s)=\mathrm{tan}^2\eta \times \varphi (\nu _a)`$.

The solar neutrino data analysis depends on four parameters: the oscillation parameters $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$, the total 8B neutrino flux $`f_{\mathrm{B},\mathrm{total}}`$ \[see eq. (16)\], and the sterile component characterized by $`\mathrm{sin}^2\eta `$. We can also characterize the sterile contribution to the total flux in terms of $`f_{\mathrm{B},\mathrm{sterile}}`$, which is defined as

$$f_{\mathrm{B},\mathrm{sterile}}\frac{\varphi (^8\mathrm{B},\mathrm{sterile})}{\varphi (^8\mathrm{B})_{\mathrm{BP00}}}.$$
(20)
We obtain from eqs. (19) and (18) the relation

$$f_{\mathrm{B},\mathrm{sterile}}=f_{\mathrm{B},\mathrm{total}}\mathrm{sin}^2\eta \left(1P_{ee}(E,\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,\mathrm{sin}^2\eta )\right).$$
(21)
Thus for a given set of parameters $`\mathrm{\Delta }m^2`$, $`\mathrm{tan}^2\theta `$, $`f_{\mathrm{B},\mathrm{total}}`$ and $`\mathrm{sin}^2\eta `$, the inferred sterile flux parameter $`f_{\mathrm{B},\mathrm{sterile}}`$ is an energy dependent, and time dependent, function. However, the energy and time dependences are mild since, for the LMA solution, $`P_{ee}`$ is nearly constant for the 8B neutrino energies that are accessible to direct measurement.

In this paper, we account for the mild energy and time dependences by defining $`f_{\mathrm{B},\mathrm{sterile}}`$ as the average sterile flux parameter inferred for the SNO experiment, i.e.,

$$f_{\mathrm{B},\mathrm{sterile}}=f_{\mathrm{B},\mathrm{total}}\mathrm{sin}^2\eta \left(1P_{ee}(E,\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,\mathrm{sin}^2\eta )_{\mathrm{SNO}}\right)$$
(22)
where $`P_{ee}(E,\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,\mathrm{sin}^2\eta )_{\mathrm{SNO}}`$ is

$$P_{ee}(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta )_{\mathrm{SNO}}=\frac{𝑑E_\nu \varphi ^{\mathrm{SSM}}({}_{}{}^{8}\mathrm{B},E_\nu )\sigma _e(E_\nu )P_{ee}(E_\nu ,\mathrm{\Delta }m^2,\mathrm{tan}^2\theta )}{R_{\mathrm{SNO}}^{\mathrm{SSM}}}.$$
(23)
Here $`P_{ee}`$ is the average survival probability and $`R_{\mathrm{SNO}}^{\mathrm{SSM}}`$ is the predicted SSM  rate for the CC SNO measurement  in the absence of neutrino oscillations,

$$R_{\mathrm{SNO}}^{\mathrm{SSM}}=𝑑E_\nu \varphi ^{\mathrm{SSM}}({}_{}{}^{8}\mathrm{B},E_\nu )\sigma _e(E_\nu ).$$
(24)
The cross section $`\sigma _e(E_\nu )`$  is the weighted average cross-section for the CC reaction, $`\nu _e+{}_{}{}^{2}He^++n+p`$, for the neutrino energy $`E_\nu `$ . The average includes the experimental energy resolution function, $`\mathrm{Res}(T,T^{})`$, where $`T(T^{})`$ is the measured (true) recoil kinetic energy of the electron. Thus

$$\sigma _e(E_\nu )=_{T_{\mathrm{min}}}^{T_{\mathrm{max}}}𝑑T_0^{T_{\mathrm{max}}^{}{}_{}{}^{}(E_\nu )}𝑑T^{}\mathrm{Res}(T,T^{})\frac{d\sigma _e(E_\nu ,T^{})}{dT^{}}.$$
(25)
The lower limit, $`T_{\mathrm{min}}`$, in the integral in eq. (25) is taken here to be the threshold, 5 MeV, used by the SNO Collaboration in ref. . The calculated value for the CC rate is not sensitive to the assumed value of $`T_{\mathrm{max}}`$, as long as $`T_{\mathrm{max}}17`$ MeV.

We construct $`\chi _{\mathrm{solar}}^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,f_{\mathrm{B},\mathrm{total}},\mathrm{sin}^2\eta )`$ in four–dimensional neutrino parameter space. The results shown in figs. 3 and 4 are the two dimensional allowed regions in $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta `$ that are obtained after marginalizing $`\chi _{\mathrm{solar}}^2`$ and $`\chi _{\mathrm{global}}^2`$ \[defined in eq. (17)\] over $`f_{\mathrm{B},\mathrm{total}}`$ and $`\mathrm{sin}^2\eta `$. The allowed domain of the oscillation parameters at a given CL is defined from the condition

$$\chi _{\mathrm{marginalized}}^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta )\chi _{min}^2+\mathrm{\Delta }\chi ^2(2,\mathrm{CL})$$
(26)
where $`\chi _{\mathrm{min}}^2=81.2`$ is the global minimum for the solar plus KamLAND analysis (for $`89=80+134`$ dof) in the four-dimensional neutrino parameter space and $`\chi _{\mathrm{marginalized}}^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta )`$ is obtained by minimizing with respect to $`f_{\mathrm{B},\mathrm{total}}`$ and $`\mathrm{sin}^2\eta `$ either $`\chi _{\mathrm{solar}}^2`$ or $`\chi _{\mathrm{global}}^2`$ for each value of $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$.

Conversely, in section 4 (and 5) we determine $`f_{\mathrm{B},\mathrm{total}}`$ (and the sterile contribution) by minimizing $`\chi _{\mathrm{global}}^2`$ with respect to $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,`$ and $`\mathrm{sin}^2\eta `$ (or $`f_{\mathrm{B},\mathrm{total}}`$). The allowed ranges for each of these parameters at a given CL are obtained from the corresponding condition for 1 dof. For instance, the allowed range of $`f_{\mathrm{B},\mathrm{total}}`$ is obtained by applying condition eq. (27), and similarly for $`\mathrm{sin}^2\eta `$ (or $`f_{\mathrm{B},\mathrm{sterile}}`$).

## 4 The total 8B solar neutrino flux

One can determine the allowed range of the total 8B neutrino flux using the results of the KamLAND reactor neutrino experiment and the results of solar neutrino experiments ref. . The conceptually simplest approach is to use just the KamLAND determination of the neutrino oscillation parameters (cf. figure 1) together with the measured CC event rate in the SNO  experiment to infer the total 8B solar neutrino flux. However, this KamLAND-only approach, while very intuitive, does not make use of the greatly reduced allowed regions obtained when solar and KamLAND data are combined (cf. figure 3 and figure 4). With the present accuracy of the KamLAND data set, the KamLAND-only determination yields an imprecise value for the total 8B neutrino flux ($`1\sigma `$ uncertainty $`30`$%).

Therefore, we present in this section the more accurate total 8B solar neutrino flux that is derived from the four–dimensional global analysis of solar plus KamLAND neutrino data (see definition of $`\chi _{\mathrm{global}}^2`$ in eq. (17) and discussion in section 3.3). Our results for $`f_{\mathrm{B},\mathrm{total}}`$ are summarized in columns three and four of table 3.

In determining $`f_{\mathrm{B},\mathrm{total}}`$, we can use all of the available solar neutrino data as well as the KamLAND measurements. This method is much more computationally demanding than the more intuitive KamLAND-only approach since in the global approach one includes, in addition to the KamLAND data, all of the 80 solar measurements and their associated (sometimes correlated) uncertainties. Formally, we calculate (as in section 3) a global $`\chi ^2`$, see the definition of $`\chi _{\mathrm{global}}^2`$ defined in eq. (17).

We obtain the allowed range of $`f_{\mathrm{B},\mathrm{total}}`$, the total 8B solar neutrino flux in units of the BP00 predicted flux \[see definition in eq. (16)\] by marginalizing $`\chi _{\mathrm{global}}^2`$ with respect to $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta `$, and $`\mathrm{sin}^2\eta `$. To be explicit, we minimize $`\chi _{\mathrm{global}}^2`$ for each value of $`f_{\mathrm{B},\mathrm{total}}`$ with respect to $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta `$, and $`\mathrm{sin}^2\eta `$. The corresponding range of $`f_{\mathrm{B},\mathrm{total}}`$ at a given CL is obtained from the condition

$$\chi _{\mathrm{marginalized}}^2(f_{\mathrm{B},\mathrm{total}})\chi _{min}^2+\mathrm{\Delta }\chi ^2(1,\mathrm{CL}).$$
(27)
We find, after extensive numerical explorations of the four-dimensional space of $`\chi _{\mathrm{global}}^2`$, that the current best empirical estimate for $`f_{\mathrm{B},\mathrm{total}}`$ is

$$f_{\mathrm{B},\mathrm{total}}=1.00(1\pm 0.06)(1\sigma ).$$
(28)
The result given in eq. (28) agrees with the standard solar model prediction .

In table 3, we show the ranges for $`f_{\mathrm{B},\mathrm{total}}`$ for within each of the LMA islands in figure 4.

## 5 The sterile component of the 8B solar neutrino flux

The sterile component of the 8B solar neutrino flux can be determined directly by subtracting the active 8B neutrino flux, which is given by the SNO NC measurement , from the total 8B neutrino flux determined in section 4 (cf. refs. ). This subtraction can be made in a model independent fashion using the oscillation parameters obtained from the analysis of the KamLAND data alone. However, this direct method, while conceptually simple, yields an imprecise determination of the sterile component of the 8B neutrino flux ($`1\sigma `$ uncertainty of order 18%).

The more powerful approach is to include all of the solar plus KamLAND measurements (cf. section 4). In this approach, we first evaluate the global $`\chi ^2`$, $`\chi _{\mathrm{global}}^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,f_{\mathrm{B},\mathrm{total}},\mathrm{sin}^2\eta )`$, which is defined in eq. (17). Not only are all of the relevant data used in calculating $`\chi _{\mathrm{global}}^2`$, but also this method properly includes the correlated uncertainties. Thus calculating $`\chi _{\mathrm{global}}^2`$ is the preferred statistical method as well as the more powerful (but computationally intensive) method.

As discussed in refs. , the allowed sterile component in the 8B neutrino flux is strongly correlated with the allowed departure from the standard solar model prediction of the total 8B neutrino flux. We illustrate this correlation in figure 6.

We show in figure 6 the allowed regions in the $`f_{\mathrm{B},\mathrm{total}}\mathrm{sin}^2\eta `$ parameter space after marginalizing with respect to $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$. For the sake of clarity, we show separately the results of marginalizing for both of the LMA islands shown in figure 4.

In figure 7, we display the corresponding allowed regions in terms of $`f_{\mathrm{B},\mathrm{total}}f_{\mathrm{B},\mathrm{sterile}}`$ (after marginalizing with respect to $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$).

Finally, marginalizing with respect to $`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta ,\mathrm{and}f_{\mathrm{B},\mathrm{total}}`$, we obtain the allowed range of the active-sterile admixture

$$\mathrm{sin}^2\eta 0.13(0.52),$$
(29)
at 1$`\sigma `$ (3$`\sigma `$).

Correspondingly for $`f_{\mathrm{B},\mathrm{sterile}}`$, we show in columns five and six of table 3 the best-fit values and the 3$`\sigma `$ ranges for the local minima of the two currently allowed LMA islands. The $`1\sigma `$ range from this global analysis is

$$f_{\mathrm{B},\mathrm{sterile}}=0.00_{0.0}^{+0.09}.$$
(30)
We use eq. (30) as our best-estimate of the sterile contribution to the observed 8B solar neutrino flux.

Is the limit on $`f_{\mathrm{B},\mathrm{sterile}}`$ given in eq. (30) (or on $`\mathrm{sin}^2\eta `$ in eq. (29) changed significantly if a more general mixing of sterile neutrinos is considered? We have investigated numerically cases in which a second $`\mathrm{\Delta }m^2`$ affects solar neutrino oscillations. We have adopted as a paradigm the model discussed in ref. , in which there is a third neutrino mass scale, $`\mathrm{\Delta }m^210^6`$eV2, related to a sterile neutrino. We find that, within the accuracy indicated by eq. (30), our best-estimate of the sterile contribution (and the associated uncertainty) of the 8B solar neutrino flux is not affected by the presence of this lower mass scale. The main effect of this low-mass sterile neutrino scenario is that the sterile neutrinos add an extra energy dependence to the survival probability and to the sterile fraction. So the limit given in eq. (30) for the 8B solar neutrino flux will not apply to the sterile component of the lower energy $`pp`$ and 7Be solar neutrino fluxes. The sterile component for the $`pp`$ and 7Be neutrino fluxes must be determined by a combination of CC and $`\nu e`$ scattering experiments at low energies ($`<1`$ MeV).

## 6 Summary

The allowed oscillation regions that are determined by including the initial KamLAND results  together with the previously known solar neutrino results (see figure 4) populate 6 orders of magnitude less area in neutrino oscillation parameter space ($`\mathrm{\Delta }m^2,\mathrm{tan}^2\theta `$) than the allowed oscillation regions that are determined with just solar neutrino data (see figure 3). The LMA solution is now the only allowed solution at $`3\sigma `$. The LOW, Vacuum, and SMA solutions are disfavored at, respectively, $`4.7\sigma `$, $`4.8\sigma `$, and $`6.1\sigma `$. All alternative explanations to the solar neutrino problem that do not also predict a deficit in KamLAND are now excluded at 3.6$`\sigma `$. Combining solar neutrino data with the KamLAND measurements, non-standard solar model explanations of the solar neutrino problem that do not include new particle physics (e.g., neutrino oscillations) are excluded at the $`7.9\sigma `$ CL. Maximal mixing, $`\mathrm{tan}^2\theta =1`$, is disfavored at $`3.1\sigma `$. The results in this paragraph are described in sec 3.1.1.

At the present level of accuracy, matter effects are not important in determining the allowed oscillation regions permitted by the KamLAND experiment. However, when the precision of the experiment is improved, matter effects must be included (see sec. 2.3). Within the currently allowed oscillation regions, the maximum change is 4% in the calculated event rate for the KamLAND experiment for a specific value of $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$.

The LMA prediction of the $`\nu e`$ scattering rate for 7Be solar neutrinos is precise, $`0.64\pm 0.03`$ within the $`1\sigma `$ (2 dof)currently allowed oscillation region (see section 3.1.2), in units of the rate predicted by the standard solar model and assuming no flavor conversion. The estimated uncertainty in the predicted standard solar model 7Be solar neutrino flux, $`\pm 10`$, is much larger than the current uncertainty in the predicted scattering rate due to neutrino oscillation parameters. Thus BOREXINO and KamLAND can provide a crucial test of solar model predictions by measuring accurately the 7Be $`\nu e`$ scattering rate.

Table 2 of section 3.1.2 summarizes the predicted rates for $`\nu e`$ scattering by $`pp`$ as well as 7Be ( and 8B) neutrinos. Within the $`1\sigma `$ (2 dof) allowed region, the predicted $`pp`$ scattering rate is $`\pm 2`$% due to imperfectly known oscillation parameters, which is about twice the currently estimated uncertainty in the solar model prediction of the $`pp`$ neutrino flux. If our current ideas regarding the accuracy of the solar model predictions are valid and if the popular neutrino oscillation framework is correct, then a measurement of the $`pp`$ solar neutrino will provide a precise determination of the mixing angle $`\theta `$.

It will be difficult with the existing SNO detector  to measure the predicted SNO CC day-night asymmetry: $`A_{\mathrm{N}\mathrm{D}}(\mathrm{SNO}\mathrm{CC})=3.3_{0.6}^{+0.7}\%(1\sigma )`$ (see eq. 9).

The event rate measured at KamLAND, (eq. 1) lies within the 1-sigma pre-KamLAND expectations based on purely solar neutrino experiments (eq. 2). This excellent agreement indicates that the main features of solar neutrinos are beginning to be well understood.

The most accurate method to determine the total 8B flux is to combine together all of the solar as well as KamLAND experimental data. Using this computationally intensive method, we find (see section 4) that the total 8B solar neutrino flux is $`1.00\left(1\pm 0.06\right)`$ in units of the standard solar model predicted flux .

The sterile component of the 8B neutrino flux can be determined similarly by using all of the available solar and KamLAND data. The allowed range of the active-sterile admixture is $`\mathrm{sin}^2\eta 0.13`$. Correspondingly, our best estimate for the sterile component of the 8B solar neutrino flux is $`0.00_{0.00}^{+0.09}`$ in units of the standard solar model predicted flux (see section 5).

Our results for the total 8B solar neutrino flux, and for the sterile component of this flux, are summarized in table 3. The present determination of the total 8B solar neutrino flux reduces the $`1\sigma `$ experimental uncertainty in this quantity by a factor of two with respect to the previous most accurate determination . The experimental bounds on the total 8B neutrino flux are now more than a factor of two smaller that the uncertainties in the predicted standard solar model flux . This is the first ocassion on which the experimental errors on a solar neutrino flux are smaller than the solar model uncertainties. The uncertainty in the sterile neutrino component of the 8B neutrino flux represents a factor of two improvement over the previous most stringent bound.

We conclude with an important methodological question. Do the results given in this paper depend in a significant way upon details of the analysis of the KamLAND data? This is a natural question to ask since the KamLAND results represent a crucial but low statistics addition to the solar neutrino data. Fortunately, different analysis strategies lead to essentially the same quantitative inferences regarding solar neutrino properties and predictions.

We have studied the robustness of our quantitative conclusions to different plausible variations in the statistical analysis of the KamLAND data. Most of these variations made no discernible difference in the final answers that were obtained from the global solar plus KamLAND analysis. For quantities such as the total 8B neutrino flux, the sterile component of the 8B neutrino flux, and the predictions for observables in future experiment, there were no significant differences. For specificity, we report here the results of a $`\chi ^2`$ analysis of the KamLAND data which led to the largest changes of any of the alternative analysis strategies we adopted. In this $`\chi ^2`$ analysis, we used the first 5 energy bins above 2.6 MeV visible energy as presented by the KamLAND collaboration but collected all the events above the fifth bin ($`E>4.725`$ MeV) into a single bin (or into two bins) in order to have enough events in each bin to assume Gaussian errors. We found that the allowed regions of the KamLAND-only analysis are slightly modified, but in ways that would be noticed only by the trained and expert eye (not to the casual reader). However, the effect of using a standard $`\chi ^2`$ with Gaussian errors is even less important for determining the global, solar plus KamLAND, allowed regions. In particular, this alternative analysis changes $`f_{\mathrm{B},\mathrm{total}}`$ by $`\pm 0.01`$ ($`1\sigma `$) and $`f_{\mathrm{B},\mathrm{sterile}}`$ by $`\pm 0.005`$ ($`1\sigma `$). The maximum change in the range of predicted solar neutrino observables in table 2 occurs for the $`R(^8\mathrm{B})`$ and amounts to $`2`$% at 1$`\sigma `$ in the largest allowed value. For all other observables, the changes are less than $`2`$% in the $`3\sigma `$ range.

We are grateful E. Akmedov and M. Smy for helpful comments and to G. Feldman for asking the important question of how large are the matter effects in the KamLAND experiment. JNB and CPG acknowledge support from NSF grant No. PHY0070928. MCG-G is supported by the European Union Marie-Curie fellowship HPMF-CT-2000-00516. This work was also supported by the Spanish DGICYT under grants PB98-0693 and by Generalitat Valenciana CTIDIB/2002/24.


format: text

text:
# A Diffusion-Based Approach to Geminate Recombination of Heme Proteins with Small Ligands

## I Introduction

The binding reactions between myoglobin ($`Mb`$) or hemoglobin ($`Hb`$) and small ligands ($`NO`$, $`O_2`$ or $`CO`$) are objects of extensive investigations during many decades because of the great functional importance of the heme proteins for living systems . In the investigations a special attention is paid to the heme-ligand recombination process going after fast photodissociative bond breaking between a ligand molecule and an ion $`Fe++`$ located in the center of a heme ($`Fe`$-protoporphyrin IX complex). The kinetic study of the postphotodissociative recombination allows to obtain detailed information on the protein-ligand interaction mechanism, the protein structure, the allosteric effect and the medium influence on the recombination efficiency (see, for example, references ).

The heme is well-wrapped in protein helixes, which prevent the iron from the solvent and hinder the ligand migration through protein matrix. On a sufficiently long time scale (at $`t100`$ ns in the case of $`Hb`$) after dissociation, when the ligand is not managed to leave the protein and to move significantly away from the parent heme, the recombination is a monomolecular reaction designated usually as a geminate recombination (GR) . Schematically, the GR can be written as :

$$A\stackrel{k_{reb}}{}B\stackrel{k_{esc}}{}\{C_1,\mathrm{},C_n\}$$
(1)
where $`A`$ is the bound heme-ligand state. The substates $`B`$ and $`\{C_1,\mathrm{},C_n\}`$ form the unbound state. Each of these substates corresponds to ligand localization in an individual cavity of protein. The substate $`B`$ answers to the residence of ligand inside the heme pocket (the cavity nearest to the iron on the distal side of heme). The rate constants $`k_{reb}`$ and $`k_{esc}`$ specify two competing processes: the irreversible heme-ligand rebinding for the ligand localized inside the heme pocket (that is, the transition from the substate $`B`$ to the state $`A`$) and the migration of unbound ligand between the heme pocket and other protein cavities (the transitions between $`B`$ and the substates $`\{C_1,\mathrm{},C_n\}`$). Immediately after photodissociation the unbound ligand is in the substate $`B`$. Therefore the quantity $`k_{esc}`$ can be associated with ligand escape from the pocket. In general, the GR is essentially determined by the specificities of heme-ligand interaction (including the spin restriction effect, the position and the orientation of ligand with respect to the heme plane)  and the effect of residues surrounding the heme . Important factors for the heme-ligand rebinding are also the state of tertiary  and quaternary  structures of protein, the conformation transitions in protein  and the solvent impact . As a consequence, the kinetic curve (that is, the probability $`P\left(t\right)`$ to find the heme in unbound state) of GR is a non-exponentially decaying function of time . After realization of the geminate stage a portion $`P_s`$ of the hemes remains in unbound state: $`P_s0.01`$ for $`NO`$, $`P_s0.1÷0.2`$ for $`O_2`$ and $`P_s0.5÷1.0`$ for $`CO`$. The quantity $`P_s`$ characterizes the efficiency of ligand escape from the protein to the solvent.

Molecular dynamics simulations  show that the movement of unbound $`NO`$, $`O_2`$ or $`CO`$ ligands in heme protein can be associated both with ligand trapping for a significant time in individual cavities and with rare jump-like transitions between adjacent cavities. It implies a fast establishment of equilibrium for the probability distribution of ligand within individual cavities. The establishment occurs on a time scale comparable to the mean time interval $`\tau _w`$ between the collisions of ligand with cavity walls. At room temperatures the time $`\tau _w`$ lies in the subpicosecond range ($`\tau _w0.1`$ ps for $`NO`$ in the heme pocket of $`Mb`$ ). The ligand redistribution between protein cavities is observed on a longer time scale ranging from several tens of picoseconds ($`40`$ ps for $`NO`$ in $`Mb`$ ) up to several tens of nanoseconds ($`50`$ ns for $`CO`$ in $`Hb`$ ). Unfortunately, in practice the detailed molecular dynamics simulation can not be implemented to the GR due to enormous computational efforts.

In the study we apply an alternative approach based on the diffusion approximation to ligand migration in protein. Such an approximation is valid for times $`t\tau _w`$ when the deterministic nature of ligand motion can be ignored. Here the interval $`\tau _w`$ can be recognized as a correlation time. The diffusion-like character of ligand migration in the heme proteins can be a reason of the non-exponential temporal dependence for the probability $`P(t)`$ . For instance, a two-dimensional diffusion is demonstrated for $`CO`$ in $`Mb`$  to explain the power-law kinetics to be observed in the experiment. Generally, reaction (1) can be represented in a three-dimensional diffusion approximation by equation

$$\frac{n}{t}=\left(Dn\right)R_{reb}n$$
(2)
with the diffusion coefficient $`D=D(x,y,z)`$. The quantity $`n=n(x,y,z,t)`$ is the probability density of unbound ligand in the protein. The stepwise function $`R_{reb}=R_{reb}(x,y,z)`$ specifies the heme-ligand rebinding and equals to $`k_{reb}`$ inside the heme pocket or to zero otherwise.

In order to solve diffusion equation (2) and to follow the evolution of GR we use a simple model proposed recently in . The model reproduces dynamics of random walk of particle in porous media (such, for instance, as glass-like matrices ) and takes into account an initial retention of ligand inside the heme pocket (that is, in the substate $`B`$). In the absence of heme-ligand rebinding the substate $`B`$ is realized at times $`t<\tau _{esc}`$ ($`\tau _{esc}=1/k_{esc}`$ is the time of ligand escape from the heme pocket to others cavities). Only on a longer time scale ($`t>\tau _{esc}`$) the ligand succeeds to leave the pocket and to migrate over protein cavities. Due to the diffusion nature of the migration the time $`\tau _{esc}`$ can be specified in terms of the diffusion coefficient $`D`$.

The approach is implemented with the help of a numerical simulation where the unbound ligand is represented by a structureless particle. For simplicity, in the simulation we make some assumptions. The ligand migration is assumed to be restricted to the distal side of heme. The ligand motion (realized on a short time scale $`t\tau _w`$) inside the heme pocket is represented by a unforced displacement of the particle within a restricted hemispheric region of space. At $`\tau _wt\tau _{esc}`$ the ligand trajectories are effectively mixed in the configurational space, resulting in a homogeneous distribution for the ligand inside the cavities. Hence, the probability of irreversible heme-ligand rebinding is accepted to be uniform for the whole heme pocket. We take into account also that on the time scale $`t\tau _w`$ the fast intracavity displacements of ligand for the substates $`\{C_1,\mathrm{},C_n\}`$ do not influence essentially on the GR kinetics and can be ignored in the simulation. Therefore the ligand displacement exterior to the heme pocket is simulated as a random walk (that is, as a Brownian-like motion) of the particle outside the hemispheric region. This walk is a spatially homogeneous diffusion with the diffusion coefficient $`D`$. We neglect also the structural transformations (such as a shift of the iron with respect to the porphyrin ring plane) at the conformational transition of protein between the unliganded and liganded states. According to the model, the temporal behavior for the probability $`P(t)`$ can be specified in terms of two parameters: the time $`\tau _{esc}`$ and the time $`\tau _{reb}=1/k_{reb}`$ of heme-ligand rebinding. The description of the model is represented in Section 2.

In order to demonstrate the usefulness of such an approach to the GR of heme proteins we apply the model to the analysis of available experimental data. We analyze the measured recombination kinetics and the efficiency for a postphotodissociative GR of human hemoglobin $`HbA`$ . These measurements were carried out at various $`pH`$ values of the solution. Here we determine the times $`\tau _{reb}`$ and $`\tau _{esc}`$ as functions of $`pH`$ and estimate the influence of solution properties on the heme-oxygen rebinding, the migration of oxygen molecule in hemoglobin and the efficiency of oxygen escape from the protein. The association of the times $`\tau _{reb}`$ and $`\tau _{esc}`$ with the time of a bimolecular recombination process for hemoglobin is analyzed. The results of simulation and their analysis are represented in Section 3.

## II Diffusion-based model of geminate recombination

The movement of unbound ligand is considered in a Cartesian coordinate system $`xyz`$ attached rigidly to the heme group of atoms. The system origin is superposed on an iron atom located in the middle of heme porphyrin ring. The $`x`$ and $`y`$ axes are aligned with the heme plane. The positive direction for the $`z`$ axis corresponds to the distal side of heme. The ligand migration in protein is simulated as a probability redistribution for the ensemble of structureless particles over a three-dimensional hemispheric space with $`z>0`$. As in , in our simulation the heme pocket is represented by a hemispheric region (designated here as a cage) of radius $`\rho `$. At an initial time instant the particle is uniformly distributed inside the cage.

The individual particle to be exposed to a sequence of $`\delta `$-shaped uncorrelated kicks executes a random walk in the space. As for the Brownian particle, each kick results in an abrupt change in the particle velocity. Between the kicks the particle is in unforced motion. On a time interval $`\mathrm{\Delta }t_k=t_{k+1}t_k`$ ($`t_k`$ is the time instant of action for $`k`$-th kick) between adjacent kicks the particle is specified by the velocity $`𝐯_k`$ and the length $`L_k`$ of free path (note that $`\mathrm{\Delta }t_k=L_k/\left|𝐯_k\right|`$ ). Then the radius vector $`𝐫(t_{k+1})`$ of particle for the time point of $`k+1`$-th kick can be obtained from iteration procedure

$$𝐫\left(t_{k+1}\right)=𝐫\left(t_k\right)+\frac{L_k𝐯_k}{\left|𝐯_k\right|}$$
(3)
where the radius vector $`𝐫(t_k)`$is given for the time instant of $`k`$-th kick. The projections $`v_{j,k}`$ ($`j=x,y,z`$) of the velocity $`𝐯_k`$ onto the coordinate axes and the length $`L_k`$ are accepted to be independent random quantities, new values of which are generated at each kick. The quantities $`v_{j,k}`$ is obtained from the Maxwell distribution

$$P_M\left(v_{j,k}\right)=\sqrt{\frac{m}{2\pi kT}}\mathrm{exp}\left(\frac{mv_{j,k}^2}{2kT}\right)$$
(4)
Here $`m`$ is the particle mass and $`T`$ is a protein temperature. At an attainment of the $`z=0`$ plane bounding the space, a new particle velocity with $`v_z>0`$ is regenerated in accordance with distribution (4).

The choice of free path length is dictated by the particle location in the space. Within the hemispheric cage the particle displacement is unforced and the particle undergoes no kicks. The length $`L_k`$ is determined then from the ballistic trajectory of particle between the cage boundaries. In this case the length is comparable to the cage size $`\rho `$. We accept here that the mean time $`\tau _h=\mathrm{\Delta }t_{h,k}`$, during which the particle crosses the cage, can be associated with the time interval $`\tau _w`$ between the collisions of ligand with heme pocket walls: $`\tau _h\tau _w`$.

Exterior to the cage, the particle is exposed to uncorrelated kicks. The absence of correlation between the kicks implies that the quantity $`L_k`$ is distributed according to the exponential law:

$$P\left(L_k\right)=\frac{1}{\lambda }\mathrm{exp}\left(\frac{L_k}{\lambda }\right)$$
(5)
where $`\lambda =L_k`$ is the mean length of free path for the particle displacement outside of the cage. The mean time $`\tau _c`$ between adjacent kicks and the length $`\lambda `$ are related to the diffusion coefficient $`D=L_k^2/6\tau _c`$ by equations:

$$\tau _c=\frac{6mD}{\pi kT}$$
(6)
$$\lambda =3D\sqrt{\frac{2m}{\pi kT}}$$
(7)
Thus, the spatial displacement of the particle is obtained from iterative equation of motion (3) and depends on the random sampling of variables $`v_{x,k}`$, $`v_{y,k}`$, $`v_{z,k}`$ and $`L_k`$, the statistical distributions for which are specified by the parameters $`m/T`$, $`D`$ and $`\rho `$. As mentioned above, under the conditions typical for the heme proteins (that is, the temperature, the ligand mass and the distinctive sizes of heme pocket) the times $`\tau _c`$ and $`\tau _h`$ to be accepted here as correlation times are negligibly short as compared to the characteristic times of GR. The length $`\lambda `$ is essentially small as against the size $`\rho `$ of hemispheric cage. Hence, the temporal behavior for the probability redistribution of ligand in heme protein can be described in terms of the diffusion-based approach.

Our model reproduces dynamics of ligand migration over protein cavities. Initially, the ligand is retained inside the heme pocket and the root-mean-square displacement $`S(t)=\sqrt{\left|𝐫(t)𝐫(0)\right|^2}`$ of ligand from the initial position does not exceed the characteristic size of the pocket. In a sense such a retention is analogous to the so called cage-effect to be observed for single atoms or small molecules in porous glass-like matrices . The time scale, on which the retention is realized, is limited by a time point $`\tau _{esc}`$. This time is a lifetime for the ligand inside the heme pocket in the absence of rebinding and specifies thereby a ligand escape from the pocket. Only on a longer time scale (when the ligand succeeds to leave the pocket and to migrate over the protein) the ligand displacement $`S(t)`$ starts to increase significantly. According to the model, we associate the time $`\tau _{esc}`$ with the time of particle localization in the hemispheric cage. In the simulation the particle displacement $`S(t)`$ does not exceed the cage radius $`\rho `$ on the short time scale $`t<\tau _{esc}`$. At longer times the quantity $`S(t)`$ increases with time. Due to the diffusion nature of the particle displacement the increase in $`S(t)^2`$ is a linear function of time and $`S(t)^26Dt`$ at $`t\tau _{esc}`$. The relation between the time $`\tau _{esc}`$ and the diffusion coefficient $`D`$ can be then determined from the requirement $`\rho ^2S(\tau _{esc})^2=6D\tau _{esc}`$:

$$\tau _{esc}=\frac{\rho ^2}{6D}$$
(8)
Fig. 1 demonstrates a typical temporal dependence for the relative particle displacement $`S(t)^2/\rho ^2`$ simulated within the framework of our model. The displacement $`S(t)`$ is shown in the figure to be constant ($`S(t)\rho `$) at $`\tau _ht\tau _{esc}`$. On a longer time scale ($`t\tau _{esc}`$) the quantity $`S(t)`$ approaches asymptotically the diffusion law: $`S(t)^26Dt=t\rho ^2/\tau _{esc}`$. Notice that for the time scale $`\tau _ht\tau _{esc}`$ the temporal behavior of relative displacements $`S(t)/\rho `$ is specified by the only parameter $`\tau _{esc}`$. In the following, we will adjust the parameter $`\tau _{esc}`$ in the simulation. For definiteness, this adjustment will be carried out by means of variation in the diffusion coefficient $`D`$. The particle mass $`m`$, the temperature $`T`$ and the cage radius $`\rho `$ will take fixed values typical for the ligand and the protein.

The heme-ligand rebinding is accepted to be an irreversible process occurring when the ligand is localized inside the heme pocket. Therefore this process is simulated as a random ’death’ for the particle within the hemispheric cage. The particle with $`\left|𝐫(t_k)\right|<\rho `$ is ’obliterated’ if $`\xi _k\mathrm{\Delta }t_{h,k}/\tau _{reb}`$. Here $`\xi _k`$ is a random quantity to be generated for each period $`\mathrm{\Delta }t_{h,k}`$ when the particle crosses the cage. The quantity $`\xi _k`$ is distributed uniformly in the interval $`[0,1]`$. The ’obliterated’ particle is excluded from the following consideration.

The probability $`P(t)`$ to find the heme in unbound state is found as the ensemble-averaged relative number of the ’non-obliterated’ particles at a time instant $`t`$. In contrast to the relative displacement $`S(t)/\rho `$, the temporal behavior for the probability $`P(t)`$ depends not only on the diffusion properties, but on the rate of heme-ligand rebinding as well. Hence, the behavior of $`P(t)`$ can be specified in terms of the times $`\tau _{esc}`$ and $`\tau _{reb}`$. In general, the probability $`P(t)`$ is a monotonously decreasing function of time, which approaches asymptotically a steady value $`P_s`$ at $`t\mathrm{}`$. This value gives a portion of the hemes remaining in unbound state after realization of GR. As in diffusion equation (2), in our model the quantity $`P_s`$ is a function dependent merely on the ratio between $`\tau _{esc}`$ and $`\tau _{reb}`$. The analysis of simulated data shows that for a wide range of values $`\tau _{esc}`$ and $`\tau _{reb}`$ satisfying the requirement $`\tau _{reb}/\tau _{esc}<20`$ (that is, under conditions typical for the GR) the best approximation of the dependence can be represented by relation

$$\frac{P_s}{1P_s}C_s\frac{\tau _{reb}}{\tau _{esc}}$$
(9)
where the coefficient $`C_s`$ is obtained from mean square fitting. At $`m=32`$ amu, $`T=300K`$ and $`\rho =4\AA `$ the fitting gives a value $`C_s=0.43`$. The inset of Fig. 1 demonstrates a good agreement between approximation (9) and the simulated data.

## III Geminate recombination of human hemoglobin with oxygen

We use the described model in order to analyze available experimental data for a postphotodissociative reoxygenation of human hemoglobin $`HbA`$. The data include the measured recombination kinetics and the efficiency of oxygen escape from the protein to the solvent for the monomolecular (geminate) and bimolecular stages of recombination reaction

$$Hb\left(O_2\right)_3+O_2Hb\left(O_2\right)_4$$
(10)
going at room temperatures after fast laser-initiated breaking of a $`FeO_2`$ bond . The kinetic measurements are carried out with a time resolution of 10 ps for a time scale $`t<1.5`$ ns at different fixed $`pH`$ values for the solution. The values of $`pH`$ fall within an interval between 6.0 and 9.4: $`pH=\{6.0,6.8,7.0,7.2,7.7,8.0,8.5,9.4\}`$.

According to the model, in the analysis of reoxygenation reaction (10) the temporal decay for the probability $`P(t)`$ is interpreted as a result of two competing processes: the heme-oxygen rebinding for the oxygen molecule localized inside the heme pocket and the diffusion-like migration of the oxygen between hemoglobin cavities. Here we determine the times $`\tau _{reb}`$ and $`\tau _{esc}`$, which specify the processes. We determine the times as functions of $`pH`$ and analyze the effect of solution properties on the processes to be considered. Due to the tetramer arrangement of hemoglobin (the $`Hb`$ molecule consists of heme containing $`\alpha `$\- and $`\beta `$-chains) the observed kinetic curve represents a reoxygenation kinetics summarized over the chains. Here we make no distinction for reaction (10) between the $`\alpha `$\- and $`\beta `$-chains and determine thereby chain-averaged times.

### A Reoxygenation kinetics for hemoglobin

The analysis of reoxygenation kinetics for the hemoglobin is based on the estimation of the times $`\tau _{reb}`$ and $`\tau _{esc}`$. The times are found with the help of a numerical simulation, the iterative procedure for which is described above (see Section 2). In the simulation the masse of walking particle is accepted to equal the mass of oxygen molecule. The temperature $`T`$ is $`300K`$. The size $`\rho `$ of hemispheric cage is $`4`$ Åthat corresponds to the time $`\tau _h<1`$ ps. The times $`\tau _{reb}`$ and $`\tau _{esc}`$ are chosen from an interval of values from 0.1 up to 5 ns. The correlation time $`\tau _c`$ and the mean length $`\lambda `$ of free path are determined by relations (6) and (7). They are negligibly small in comparison with $`\tau _{reb}`$, $`\tau _{esc}`$ or $`\rho `$. The simulated dependences for the probability $`P(t)`$ are obtained from ensemble averaging for more than $`10^6`$ particles.

In the simulation the parameters $`\tau _{reb}`$ and $`\tau _{esc}`$ are so adjusted that the ensemble-averaged temporal dependence of simulated probability $`P(t)`$ is the best agreement with a measured kinetic curve. The agreement is specified by the relative root-mean-square deviation $`R`$ between the simulated and experimental curves. The simulated dependence for $`P(t)`$ is shown in Fig. 2 to reproduce well kinetic measurements on the considered time scale. The minimal deviation $`R`$ achieved in our calculations for each of the fixed $`pH`$ values does not exceed the measurement error ($`R0.01`$). Such an agreement testifies that the non-exponential dependence for $`P(t)`$ with time can be explained by a diffusion-like migration of ligand over protein matrix. Hence, the parameters $`\tau _{reb}`$ and $`\tau _{esc}`$ can be used for the analysis of the processes, which are responsible for the GR.

The influence of solution properties on the migration and the rebinding of oxygen molecule in hemoglobin is assessed from a $`pH`$ dependence for the obtained times $`\tau _{reb}`$ and $`\tau _{esc}`$. Our simulation demonstrates a significant variation in the rate of heme-oxygen rebinding with $`pH`$ (see Fig. 3). The increase of quantity $`pH`$ from 6.0 to 8.5 results in a shortening for the time $`\tau _{reb}`$ by a factor of 1.4 (from 0.5 down to 0.36 ns). With the following rise of $`pH`$ to 9.4 the parameter $`\tau _{reb}`$ appears to increase up to 0.4 ns. The minimum magnitude of $`\tau _{reb}`$ is observed at $`pH=8.5`$. Despite the considerable $`pH`$ effect for the heme-oxygen rebinding, the variation in $`pH`$ influences slightly on the oxygen escape from the heme pocket. The time $`\tau _{esc}`$ is shown in Fig. 3 to be within a range of values from 0.82 to 0.92 ns and to be weakly dependent on $`pH`$. The $`pH`$-averaged magnitude of $`\tau _{esc}`$ is approximately 0.88 ns. Notice that this magnitude is larger than $`\tau _{reb}`$ by a factor of $`2÷3`$.

The obtained values for the times $`\tau _{reb}`$ and $`\tau _{esc}`$ are in good agreement with the experimental study of the alkaline Bohr effect (the variation of the recombination rate $`1/\tau _s`$ with $`pH`$ for the bimolecular stage of GR) . The behavior of $`pH`$-dependence for the time $`\tau _{reb}`$ is demonstrated in Fig. 3 to be similar to one for the time $`\tau _s`$ of bimolecular rebinding. Such an agreement testifies that for the monomolecular GR the variation of the rebinding rate with $`pH`$ can be associated with the same structural transformation as for the bimolecular stage of reaction (10). Histidine imidazoles of $`C`$-terminal sites and $`\alpha `$-amides of $`N`$-terminal sites seem to be the aminoacid residues, which are responsible for this transformation . Specifically, in the alkaline Bohr effect the interaction between the solvent and the $`\beta 146His`$ residue (a $`C`$-terminal histidine of $`\beta `$-chain) is one of the most probable reasons for the heme structure modification and the rearrangement of neighboring aminoacid residues . Our simulation confirms that the variation in $`pH`$ can result in essential structural transformations in immediate proximity from the iron atom. The strong $`pH`$-dependence for the times $`\tau _{reb}`$ and $`\tau _s`$ is a consequence of the transformations.

The ligand penetration from the solvent into the heme pocket is shown for $`O_2`$ or $`NO`$ in $`Hb`$  to be a process restraining the rate of bimolecular recombination (10). Therefore, the similarity between the $`pH`$-dependences for $`\tau _{reb}`$ and $`\tau _s`$ testifies that the change in $`pH`$ has a slight effect on the oxygen migration in hemoglobin at the mono- and bimolecular stages of recombination (10). The weak $`pH`$-dependence for the obtained times of oxygen escape from the pocket confirms this assumption. Such a $`pH`$-invariant behavior for the oxygen migration can be interpreted by the independence of mobility for the hemoglobin side chains (which seem to be responsible for the ligand transitions between cavities of heme protein ) on $`pH`$ of the solution.

### B Efficiency of oxygen escaping from hemoglobin

The obtained times $`\tau _{reb}`$ and $`\tau _{esc}`$ are used then in order to estimate the efficiency of oxygen escape from hemoglobin as a function of $`pH`$. The efficiency is proportional to the quantum yield of photodissociation and can be associated with the portion $`P_s`$ of the hemes remaining in unbound state after realization of the geminate reoxygenation stage . In our simulation the ratio $`\tau _{reb}/\tau _{esc}`$ falls within a range of values from 0.4 up to 0.6. It implies that the portion $`P_s`$ can be determined from approximation(9). The behavior of $`pH`$-dependence for the obtained quantity $`P_s`$ agrees well with one for the measured quantum yield of photodissociation . The quantity $`P_s`$ is shown in Fig. 4 to be proportional to the apparent quantum yield $`\gamma `$ for the whole investigated scale of $`pH`$.

Notice that for the studied range of $`pH`$ values the quantity $`C_s\tau _{reb}/\tau _{esc}`$ is considerably low in comparison with 1 ($`C_s\tau _{reb}/\tau _{esc}0.15÷0.25`$) and the time $`\tau _{esc}`$ is practically constant. Therefore the relation of the portion $`P_s`$ with the time $`\tau _{reb}`$ is close to the linear law:

$$P_s=\frac{C_s\tau _{reb}}{\left(\tau _{esc}+C_s\tau _{reb}\right)}C_s\frac{\tau _{reb}}{\tau _{esc}}\tau _{reb}.$$
(11)
In our study (see Fig. 3 and Fig. 4) the $`pH`$-dependences obtained for $`P_s`$ and $`\tau _{reb}`$ are similar that testifies again that the transport properties for the oxygen molecule in hemoglobin do not depend on $`pH`$.

### C Diffusion properties of oxygen migration in hemoglobin

The analysis of X-ray diffraction data  for oxygenated and deoxygenated species of human hemoglobin (PDB ID $`1HHO`$ and $`2HHB`$, correspondingly) shows that the cage radius $`\rho `$ to be associated with the heme pocket size is a quantity ranging from 1 up to 5 Å(taking into account the Van der Waals radiuses). Hence, the diffusion coefficient $`D`$ for the oxygen migration in hemoglobin can be estimated from relation (8): $`D=\rho ^2/6\tau _{esc}0.2÷510^{11}`$ m<sup>2</sup>/s. This coefficient is intermediate to diffusion coefficients for small molecules in water ( 10<sup>-9</sup> m<sup>2</sup>/s) and solids (10<sup>-18</sup> m<sup>2</sup>/s at $`T<400K`$)

According to the diffusion law, at a time instant $`t_m`$ the root-mean-square displacement $`\sqrt{\left|𝐫(𝐭)\right|^2}`$ of the ligand from the iron is approximately equal to $`\rho \sqrt{t_m/\tau _{esc}}`$. It implies that on the completion of kinetic measurements ($`t_m=1.5`$ ns ) the oxygen remains inside the protein and is localized in immediate proximity from the heme pocket: $`\sqrt{\left|𝐫(𝐭)\right|^2}1.3\rho <7`$ Å. This conclusion is consistent with results of spectroscopy investigation of motional dynamics for $`CO`$ in $`Hb`$ .

## IV Conclusion

We have represented a simple model of the geminate recombination of heme proteins with small ligands. The model takes into account dynamic properties of ligand displacement in protein matrix. In the model the recombination is due both to the heme-ligand rebinding and to the diffusion-like migration of ligand between protein cavities. The temporal behavior for the probability $`P(t)`$ to find the heme in unbound state is specified in terms of two parameters. They are the time $`\tau _{reb}`$ of heme-ligand rebinding for the ligand inside the heme pocket and the time $`\tau _{esc}`$ of ligand escape from the pocket.

We have applied our model in order to analyze a postphotodissociative geminate reoxygenation of human hemoglobin at various $`pH`$ values of the solution. The measured kinetic curves and the efficiency of oxygen escape from the hemoglobin are well reproduced in our simulation. It testifies that the non-exponential behavior for the probability $`P(t)`$ can be explained by a diffusion-like migration of ligand over protein cavities. This conclusion is consistent with recent kinetic measurements . We believe that the theory-experiment agreement may be considered as an additional validation for the glass-like model of proteins.

Our study demonstrates also that the variation in $`pH`$ can result in considerable changes for the rate of heme-ligand rebinding. At the time, the oxygen migration in hemoglobin depends slightly on $`pH`$. We have interpreted this effect as a result of essential structural transformations in immediate proximity from the iron atom. Certainly, this conclusion demands a more detailed and thorough examination. In any case we suppose that the $`pH`$-induced modification of the initial stage of GR (if the modifications are observed) can be explained by a change in the rate of heme-ligand rebinding.