SEARCHES FOR LARGE-SCALE ANISOTROPY IN THE ARRIVAL DIRECTIONS OF COSMIC RAYS DETECTED ABOVE ENERGY OF 10¹⁹ eV AT THE PIERRE AUGER OBSERVATORY AND THE TELESCOPE ARRAY

The large-scale distribution of arrival directions of cosmic rays is an important observable in attempts to understand their origin. This is because this observable is closely connected to both their source distribution and propagation. Due to scattering in magnetic fields, the anisotropy imprinted upon the distribution of arrival directions is mainly expected at large scales up to the highest energies. Large-scale patterns with anisotropy contrast at the level of 10⁻⁴ to 10⁻³ have been reported by several experiments for energies below 10¹⁵ eV where the high flux of cosmic rays allows the collection of a large number of events (Amenomori et al. 2005; Guillian et al. 2007; Aglietta et al. 2009; Abdo et al. 2009; Abbasi et al. 2010, 2011, 2012; Aartsen et al. 2013). For energies above a few 10¹⁵ eV, the decrease of the flux with energy makes it challenging to collect the necessary statistics required to reveal amplitudes down to 10⁻³ or 10⁻². Upper limits at the level of a few percent have been obtained at ≃ 10¹⁶ eV (Curcio et al. 2013) and ≃ 10¹⁸ eV (The Pierre Auger Collaboration 2011a, 2012).

The anisotropy of any angular distribution on the sphere is encoded in the corresponding set of spherical harmonic moments a_ℓm. Although not predictable in a quantitative way at present, large-scale anisotropies might be expected from various mechanisms of propagation of cosmic rays. A non-zero dipole moment is naturally expected from propagation models leading to a cosmic ray density gradient embedding the observer. Even in the absence of a density gradient, a measurable dipole moment might result from the motion of the Earth or of the massive objects in the neighborhood of the Milky-Way relative to a possibly stationary cosmic-ray rest frame (Compton & Getting 1935; Kachelriess & Serpico 2006; Harari et al. 2010). On the other hand, excesses along a plane, for instance the super-Galactic one, would be detectable as a prominent quadrupole. The dipole and the quadrupole moments are thus of special interest, but an access to the full set of multipoles is relevant to characterize departures from isotropy at all scales. However, since cosmic ray observatories at ground level have only a partial-sky coverage, the recovering of these multipoles turns out to be nearly impossible without explicit assumptions on the shape of the underlying angular distribution (Sommers 2001). Indeed, for an angular distribution described by a multipolar expansion bounded to some degree ℓ, the multipole coefficients can only be estimated within a resolution degraded exponentially with ℓ (Billoir & Deligny 2008). In most cases, given the available statistics, only the dipole and the quadrupole moments can be estimated with a relevant resolution under the assumption that the flux of cosmic rays is purely dipolar or purely dipolar and quadrupolar, respectively. Evading such hypotheses and thus measuring the multipoles to any order in an unambiguous way requires full-sky coverage. At present, full-sky coverage can only be achieved through the meta-analysis of data recorded at observatories located in both hemispheres.

The Telescope Array and the Pierre Auger Observatory are the two largest experiments ever built to study ultra-high energy cosmic rays in each hemisphere. The aim of the joint analysis reported in this article is to search for anisotropy with full-sky coverage by combining the data of the two experiments. The data sets used for this search are described in Section 2 together with some properties and performances of the experiments relevant for this study. Special emphasis is given to the control of the event counting rate with time and local angles as well as to the respective exposures to each direction of the sky. To facilitate this first joint analysis, the energy threshold used in this report, 10¹⁹ eV, is chosen to guarantee that both observatories operate with full detection efficiency for any of the events selected in each data set. Above this energy, the respective exposure functions follow purely geometric expectations.

The analysis method to estimate the spherical harmonic moments a_ℓm is presented in Section 3, together with its statistical performance. The main challenge in combining the data sets is to account adequately for the relative exposures of both experiments. The empirical approach adopted here is shown in Section 4 to meet the challenge. Results of the estimated multipole coefficients are then presented and illustrated in several ways in Section 5 with, in particular, reports for the first time of a significance full-sky map of the overdensities and underdensities and of the angular power spectrum above 10¹⁹ eV. Several cross-checks against systematic effects are presented in Section 6, showing the robustness of the analyses. Finally, the results are discussed in Section 7, together with some prospects for future joint analyses.

2.1. The Pierre Auger Observatory

2.2. The Telescope Array

2.3. Control of the Event Rate

Control of the event rate is of critical importance. The magnitude of the spurious pattern imprinted in the arrival directions by any effect of experimental origin must be kept under control. This is essential if this magnitude is larger than the fluctuations on the anisotropy parameters intrinsic to the available statistics.

The instantaneous exposure of each experiment is not constant in time due to the construction phase of the observatories and to unavoidable dead times of detectors (e.g., failures of electronics, power supply, communication system, etc.). This translates into modulations of the event rates even for an isotropic flux. However, these dead times concern only a few detectors and are randomly distributed over operation time, so that once averaged over several years of data taking, the relative modulations of both exposure functions in local sidereal time turn out to be not larger than a few per thousand. The impact for anisotropy searches is thus expected to be negligible given the small statistics available.

In terms of local angles, the event counting rate is controlled by the relationship between the measured shower size and zenith angle used to estimate the energy by accounting for the attenuation of the showers in the atmosphere. For an isotropic flux and full efficiency, the distribution in zenith angles dN/dθ of a surface detector array is expected to be proportional to sin θ cos θ for solid angle and geometry reasons, so that the distribution in dN/dsin ²θ is expected to be uniform. Note that this distribution is quasi-invariant to large-scale anisotropies (The Pierre Auger Collaboration 2012), so that requiring the distributions dN/dsin ²θ to be uniform constitutes a well-suited tool to control the event counting rate. Both distributions are shown in Figure 1 to be indeed uniform above 10¹⁹ eV within statistical uncertainties.

Zoom In Zoom Out Reset image size

Due to the steepness of the energy spectrum, the event counting rate can also be largely distorted by systematic changes of the energy estimate with time and/or local angles. The two main effects acting in this way, namely, the atmospheric and geomagnetic effects, are by default accounted for in large-scale anisotropy searches reported by the Auger collaboration (The Pierre Auger Collaboration 2008, 2011b). This is necessary given the available statistics around 10¹⁸ eV. However, above 10¹⁹ eV, the impact of these effects is marginal given the reduced statistics. The fact that the same treatment is not implemented in the data set provided by the Telescope Array has no impact on the accuracy of the anisotropy measurements presented in this report, as is shown below.

Changes of atmospheric conditions are known to modulate the event rate as a function of time. This is because the development of an extensive air shower is sensitive to the atmospheric pressure and air density in a way which influences the measurement of the shower size at a fixed distance from the core, and consequently the measurement of the energy. To avoid the undesired variations of the event rate induced by these effects, the observed shower size has to be related to the one that would have been measured at some fixed reference values of pressure and density. Such a procedure is implemented to produce the data set recorded at the Pierre Auger Observatory (The Pierre Auger Collaboration 2008).

Assuming there is no variation of the cosmic-ray flux over a timescale of a few years, the fact that the time windows of the data sets provided by the two observatories are different has no impact on the anisotropy searches presented in this report, given that more than three years of data are considered in each set separately. With at least three years of data taking, indeed, the event rate variation at the solar timescale decouples from the one at the sidereal timescale, so that any significant modulation of the event rate of experimental origin, primarily visible at the solar timescale, would not significantly impact the event rate at the sidereal timescale (Billoir & Letessier-Selvon 2008).

Although there is no shower size correction for weather effects in the data set of the Telescope Array, it is worth noting, however, that the natural timescale at which the modulation of the event rate operates is the solar one. This means that over whole years of data taking as in the present analysis, such a modulation is expected to be partially compensated at the sidereal timescale, as just emphasized. The size of the residual effect, together with the level of the sideband effect induced by a seasonal modulation of the daily counting rate, can be checked empirically by evaluating the amplitude of the first harmonics for fictitious right ascension angles calculated by dilating the time of the events in such a way that a solar day lasts about four minutes longer (Farley & Storey 1954). The corresponding timescale is called the anti-sidereal one. A first harmonic amplitude standing out from the background noise at such an unphysical timescale would be indicative of important spurious effects of instrumental origin in the measurement of the first harmonic coefficients at the sidereal timescale. The measured values are shown in Figure 2 for each experiment, together with the respective Rayleigh distributions expected from statistical fluctuations. The amplitudes are seen to be compatible with that expected from the Rayleigh distributions. This provides support that the counting rate is not affected by spurious modulations of instrumental origin at the Pierre Auger Observatory—as expected from the corrections of the signal sizes—and at the Telescope Array as well.

Zoom In Zoom Out Reset image size

Moreover, since we aim at characterizing the arrival directions in both right ascension and declination angles, it is also important to control the event rate in terms of local angles. The geomagnetic field turns out to influence shower developments and shower size estimations at a fixed energy due to the broadening of the spatial distribution of particles in the direction of the Lorentz force. The strength of the resulting modulation of the event rate depends on the angle between the incoming direction of each event and the direction of the transverse component of the geomagnetic field. The event rate is thus distorted independently of time as a function of local zenith and azimuth angles, and thus as a function of the declination. To eliminate these distortions, the data set recorded at the Pierre Auger Observatory is produced by relating the shower size of each event to the one that would have been observed in the absence of the geomagnetic field (The Pierre Auger Collaboration 2011b). However, the impact of these corrections is shown in Section 6 to be marginal given the limited statistics above 10¹⁹ eV. That the same kind of corrections are not carried out on the data set provided by the Telescope Array is thus unimportant for the present study. This is further reinforced by the fact that geomagnetic effects are expected to be more important for a water-Cherenkov detector array which is more sensitive to muonic signal than for a scintillator one as measurements are made up to larger zenith angles.

2.4. Directional Exposures

The directional exposure ω(n) provides the effective time-integrated collecting area for a flux from each direction of the sky. The small variations of the exposure in sidereal time translate into small variations of the directional exposure in right ascension. However, given the small size of these variations, the relative modulations of the respective directional exposure functions in right ascension turn out to be less than a few 10⁻³. Given that the limited statistics currently available above 10¹⁹ eV cannot allow an estimation of each a_ℓm coefficient with a precision better than a few percent, the non-uniformities of both ω_TA and ω_Auger in right ascension can be neglected. Both functions are consequently considered to depend only on the declination hereafter.

Since the energy threshold of 10¹⁹ eV guarantees that both experiments are fully efficient in their respective zenith range [0, θ_max], the directional exposure relies only on geometrical acceptance terms. The dependence on declination can then be obtained in an analytical way (Sommers 2001) as

$\begin{equation} \omega _i(\mathbf {n})=A_i(\cos \lambda _i\,\cos \delta \,\sin \alpha _m+\alpha _m\sin \lambda _i\,\sin \delta), \end{equation} \tag{ 1 }$

where λ_i is the latitude of the considered experiment, the parameter α_m is given by

$\begin{equation} \alpha _m=\left\lbrace \begin{array}{@{}l@{\quad }l@{}} 0 ;&\xi >1,\\ \pi ;&\xi <-1,\\ \arccos \xi ;& {\rm otherwise,} \end{array}\right. \end{equation} \tag{ 2 }$

with ξ ≡ (cos θ_max − sin λ_i sin δ)/cos λ_i cos δ, and the normalization factors A_i chosen such that the integration of each ω_i function over 4π matches the (total) exposure of the corresponding experiment. The directional exposure functions ω_i(δ) of each experiment are shown in Figure 3. Given the respective latitudes of both observatories and with the maximum zenith angle used here, overall, it is clearly seen that full-sky coverage is indeed achieved when summing both functions. Also, and it will be important in the following, it is interesting to note that a common band of declination, namely, −15° ⩽ δ ⩽ 25°, is covered by both experiments.

Zoom In Zoom Out Reset image size

In principle, the combined directional exposure of the two experiments should be simply the sum of the individual ones. However, individual exposures have to be re-weighted here by some empirical factor b due to the unavoidable uncertainty in the relative exposures of the experiments:

$\begin{equation} \omega (\mathbf {n};b)=\omega _{\rm TA}(\mathbf {n})+b\,\omega _{\rm Auger}(\mathbf {n}). \end{equation} \tag{ 3 }$

Written in this way, b is a dimensionless parameter of order unity. In practice, only an estimation $\bar{b}$ of the factor b can be obtained, so that only an estimation of the directional exposure $\bar{\omega }(\mathbf {n})\equiv \omega (\mathbf {n};\bar{b})$ can be achieved through Equation (3). The procedure used for obtaining $\bar{b}$ from the joint data set is described in Section 4. In addition, although the techniques for assigning energies to events are nearly the same, there are differences as to how the primary energies are derived at the Telescope Array and the Pierre Auger Observatory. Currently, systematic uncertainties in the energy scale of both experiments amount to about 21% and 14%, respectively (Abu-Zayyad et al. 2013b; Verzi et al. 2013). This encompasses the adopted fluorescence yield, the uncertainties for the absolute calibration of the fluorescence telescopes, the influence of the atmosphere transmission used in the reconstruction, the uncertainties in the shower reconstruction, and the uncertainties in the correction factor for the missing energy. Uncovering and understanding the sources of systematic uncertainties in the relative energy scale is beyond the scope of this report and will be addressed elsewhere. However, such a potential shift in energy leads to different counting rates above some fixed energy threshold, which induces fake anisotropies. Formally, these fake anisotropies are similar to the ones resulting from a shift in the relative exposures of the experiments.¹³¹ The parameter b can thus be viewed as an effective correction which absorbs any kind of systematic uncertainties in the relative exposures, whatever the sources of these uncertainties. This empirical factor is arbitrarily chosen to re-weight the directional exposure of the Pierre Auger Observatory relatively to the one of the Telescope Array.

The observed angular distribution of cosmic rays, dN/dΩ, can be naturally modeled as the sum of Dirac functions on the surface of the unit sphere whose arguments are the arrival directions {n₁, ..., n_N} of the events,

$\begin{equation} \frac{{d}N(\mathbf {n})}{{d}\Omega }=\sum _{i=1}^N\delta (\mathbf {n},\mathbf {n}_i). \end{equation} \tag{ 4 }$

Throughout this section, arrival directions are expressed in the equatorial coordinate system (declination δ and right ascension α) since this is the most natural one tied to the Earth in describing the directional exposure of any experiment. The random sample {n₁, ..., n_N} results from a Poisson process whose average is the flux of cosmic rays Φ(n) coupled to the directional exposure ω(n) of the considered experiment,

$\begin{equation} \left\langle \frac{{d}N(\mathbf {n})}{{d}\Omega }\right\rangle =\omega (\mathbf {n})\,\Phi (\mathbf {n}). \end{equation} \tag{ 5 }$

As for any angular distribution on the unit sphere, the flux of cosmic ray Φ(n) can be decomposed in terms of a multipolar expansion onto the spherical harmonics Y_ℓm(n),

$\begin{equation} \Phi (\mathbf {n})=\sum _{\ell \ge 0}\sum _{m=-\ell }^\ell a_{\ell m}Y_{\ell m}(\mathbf {n}). \end{equation} \tag{ 6 }$

Any anisotropy fingerprint is encoded in the multipoles a_ℓm. Non-zero amplitudes in the ℓ modes contribute in variations of the flux on an angular scale ≃ 1/ℓ radians. The rest of this section is dedicated to the definition of an estimator $\bar{a}_{\ell m}$ of the multipolar coefficients and to the derivation of the statistical properties of this estimator.

With full-sky but non-uniform coverage, the customary recipe for decoupling directional exposure effects from anisotropy ones consists in weighting the observed angular distribution by the inverse of the relative directional exposure function (Sommers 2001)

$\begin{equation} \frac{{d}\tilde{N}(\mathbf {n})}{{d}\Omega }=\frac{1}{\bar{\omega }_r(\mathbf {n})}\frac{{d}N(\mathbf {n})}{{d}\Omega }. \end{equation} \tag{ 7 }$

The relative directional exposure $\bar{\omega }_r$ is a dimensionless function normalized to unity at its maximum. When the function ω (or ω_r) is known from a single experiment, the averaged angular distribution $\langle {d}\tilde{N}/{d}\Omega \rangle$ from Equation (5) is identified with the flux of cosmic rays Φ(n) times the total exposure of the experiment. In turn, when combining the exposure of the two experiments, the relationship between $\langle {d}\tilde{N}/{d}\Omega \rangle$ and Φ(n) is no longer so straightforward due to the finite resolution in estimating the parameter b introduced in Equation (3). To first order, it can be expressed as

$\begin{equation} \left\langle \frac{{d}\tilde{N}(\mathbf {n})}{{d}\Omega }\right\rangle \simeq \left\langle \frac{1}{\bar{\omega }_r(\mathbf {n})}\right\rangle \omega (\mathbf {n})\,\Phi (\mathbf {n}). \end{equation} \tag{ 8 }$

For an unbiased estimator of b with a resolution¹³² not larger than ≃ 10%, the relative differences between ${\langle }1/\bar{\omega }_r(\mathbf {n})\rangle$ and 1/ω_r(n) are actually not larger than 10⁻³ in such a way that $\langle {d}\tilde{N}/{d}\Omega \rangle$ can still be identified to Φ(n) times the total exposure to a very good approximation. Consequently, the $\bar{a}_{\ell m}$ coefficients recovered, defined as

$\begin{equation} \bar{a}_{\ell m}=\int _{4\pi }{d}\Omega \frac{{d}\tilde{N}(\mathbf {n})}{{d}\Omega }Y_{\ell m}(\mathbf {n})=\sum _{i=1}^N\frac{Y_{\ell m}(\mathbf {n}_i)}{\bar{\omega }_r(\mathbf {n}_i)}, \end{equation} \tag{ 9 }$

provide unbiased estimators of the underlying a_ℓm multipoles since the relationship $\langle \bar{a}_{\ell m}\rangle =a_{\ell m}$ can be established by propagating Equation (8) into $\langle \bar{a}_{\ell m}\rangle$.

Using the estimators defined in Equation (9), the expected resolution σ_ℓm on each a_ℓm multipole can be inferred from the second moment of ${d}\tilde{N}/{d}\Omega$ accordingly to Poisson statistics,

$\begin{eqnarray} &&\left\langle \frac{{d}\tilde{N}(\mathbf {n})}{{d}\Omega }\frac{{d}\tilde{N}(\mathbf {n}^{\prime })}{{d}\Omega ^{\prime }}\right\rangle =\left\langle \frac{1}{\bar{\omega }_r(\mathbf {n})\bar{\omega }_r(\mathbf {n}^{\prime })}\right\rangle\nonumber\\ &&\qquad\times [\omega (\mathbf {n})\omega (\mathbf {n}^{\prime })\,\Phi (\mathbf {n})\,\Phi (\mathbf {n}^{\prime })+\omega (\mathbf {n})\,\Phi (\mathbf {n})\,\delta (\mathbf {n},\mathbf {n}^{\prime })].\nonumber\\ && \end{eqnarray} \tag{ 10 }$

Once propagated into the covariance matrix of the estimated $\bar{a}_{\ell m}$ coefficients, Equation (10) allows the determination of σ_ℓm in the case of relatively small {a_ℓm}_{ℓ ⩾ 1} coefficients compared to a₀₀. Using as above the fact that ${\langle }1/\bar{\omega }_r(\mathbf {n})\rangle$ can be accurately replaced by 1/ω_r(n), the resolution parameters σ_ℓm read

$\begin{eqnarray} \sigma _{\ell m}&\simeq& \left[\frac{a_{00}}{\sqrt{4\pi }}\int _{4\pi }{d}\Omega \left\langle \frac{1}{\bar{\omega }_r^2(\mathbf {n})}\right\rangle \omega (\mathbf {n})\,Y_{\ell m}^2(\mathbf {n})\right.\nonumber \\ &&\left.+ \frac{a_{00}^2}{4\pi }\int _{4\pi }{d}\Omega \,{d}\Omega ^{\prime } \left[\left\langle \frac{1}{\bar{\omega }_r(\mathbf {n})\,\bar{\omega }_r(\mathbf {n}^{\prime })}\right\rangle\right.\right.\nonumber\\ &&\left.\left.\times\, \omega (\mathbf {n})\,\omega (\mathbf {n}^{\prime })-1\right]Y_{\ell m}(\mathbf {n})\,Y_{\ell m}(\mathbf {n}^{\prime })\right]^{1/2}. \end{eqnarray} \tag{ 11 }$

If b was known with perfect accuracy, the second term in Equation (11) would vanish, and the resolution of the a_ℓm coefficients would be driven by Poisson fluctuations only as in the case of a single experiment. However, having at one's disposal an estimation of b only, the second term reflects the effect of the uncertainty in the relative exposures of the two experiments. For a directional exposure independent of the right ascension, the azimuthal dependences of the spherical harmonics can be factorized from the whole solid angle integrations so that the whole term is non-zero only for m = 0. Its influence is illustrated in Figure 4, where the ratio between the total expression of σ_ℓ0 and the partial one, ignoring this second term inside the square root, is plotted as a function of the multipole ℓ for different resolution values on b. While this ratio amounts to ≃ 1.5 for ℓ = 1 and σ(b)/b = 3.5%, it falls to ≃ 1.1 for ℓ = 2 and then tends to 1 for higher multipole values. Consequently, in accordance with naive expectations, the uncertainty in the b factor mainly impacts the resolution in the dipole coefficient a₁₀ while it has a small influence on the quadrupole coefficient a₂₀ and a marginal one on higher order moments {a_ℓ0}_{ℓ ⩾ 3}.

Zoom In Zoom Out Reset image size

A band of declinations around the equatorial plane is exposed to the fields of view of both experiments, namely, for declinations between −15° and 25°. This overlapping region can be used for designing an empirical procedure to obtain a relevant estimate of the parameter b. The basic starting point is the following. For an isotropic flux, the fluxes measured independently by both experiments in the common band would have to be identical. The commonly covered declination band could thus be used for cross-calibrating empirically the fluxes of the experiments and for delivering an overall unbiased estimate of the parameter b. Since the shapes of the exposure functions are not identical in the overlapping region (see Figure 3), the observed fluxes are not expected to be identical in case of anisotropies. For small anisotropies, however, this guiding idea can nevertheless be implemented in an iterative algorithm, finally delivering estimates of the parameter b and of the multipole coefficients at the same time.

Let us consider a joint data set with all events detected in excess of some energy threshold. The way the individual energy scales are chosen to select a common threshold, whether using nominal energies or any cross-calibration procedure, does not matter at this stage. For anisotropies which do not vary suddenly with energy, only a reasonable starting point is required to guarantee that the anisotropy search pertains to events with energies in excess of roughly the same energy threshold for both experiments. Then, considering as a first approximation the flux Φ(n) as isotropic, the overlapping region ΔΩ can be utilized to derive a first estimate $\bar{b}^{(0)}$ of the b factor by forcing the fluxes of both experiments to be identical in this particular region. This can be easily achieved in practice, by taking the ratio of the ΔN_TA and ΔN_Auger events observed in the overlapping region ΔΩ weighted by the ratio of nominal exposures,

$\begin{equation} \bar{b}^{(0)}=\frac{\Delta N_{\rm Auger}}{\Delta N_{\rm TA}}\,\frac{\int _{\Delta \Omega }{d}\Omega \,\omega _{\rm TA}(\mathbf {n})}{\int _{\Delta \Omega }{d}\Omega \,\omega _{\rm Auger}(\mathbf {n})}. \end{equation} \tag{ 12 }$

Then, inserting $\bar{b}^{(0)}$ into $\bar{\omega }$, "zero-order" $\bar{a}_{\ell m}^{(0)}$ coefficients can be obtained. This set of coefficients is only a rough estimation, due to the limiting assumption on the flux (isotropy).

On the other hand, the expected number of events in the common band for each observatory, $\Delta n_{\rm TA}^{\rm exp}$ and $\Delta n_{\rm Auger}^{\rm exp}$, can be expressed from the underlying flux Φ(n) and the true value of b as

$\begin{eqnarray} \Delta n_{\rm TA}^{\rm exp}&=\int _{\Delta \Omega }{d}\Omega \,\Phi (\mathbf {n})\,\omega _{\rm TA}(\mathbf {n}) \nonumber \\ \Delta n_{\rm Auger}^{\rm exp}&=b\int _{\Delta \Omega }{d}\Omega \,\Phi (\mathbf {n})\,\omega _{\rm Auger}(\mathbf {n}). \end{eqnarray} \tag{ 13 }$

From Equations (13), and from the set of $\bar{a}_{\ell m}^{(0)}$ coefficients, an iterative procedure estimating at the same time b and the set of a_ℓm coefficients can be constructed as

$\begin{equation} \bar{b}^{(k+1)}=\frac{\Delta N_{\rm Auger}}{\Delta N_{\rm TA}}\,\frac{\int _{\Delta \Omega }{d}\Omega \,\bar{\Phi }^{(k)}(\mathbf {n})\,\omega _{\rm TA}(\mathbf {n})}{\int _{\Delta \Omega }{d}\Omega \,\bar{\Phi }^{(k)}(\mathbf {n})\,\omega _{\rm Auger}(\mathbf {n})}, \end{equation} \tag{ 14 }$

where ΔN_TA and ΔN_Auger as derived in the first step are used to estimate $\Delta n_{\rm TA}^{\rm exp}$ and $\Delta n_{\rm Auger}^{\rm exp}$, respectively, and $\bar{\Phi }^{(k)}$ is the flux estimated with the set of $\bar{a}_{\ell m}^{(k)}$ coefficients.

Whether this iterative procedure finally delivers unbiased estimations of the set of a_ℓm coefficients with a resolution obeying Equation (11) can be tested by Monte-Carlo simulations. 10,000 mock samples are used here, with a number of events similar to the one of the actual joint data set and with ingredients corresponding to the actual figures in terms of total and directional exposures. Under these realistic conditions, the resolution obtained on the b parameter is found to be ≃ 3.9%. The distributions of the reconstructed low-order $\bar{a}_{10}$ and $\bar{a}_{20}$ multipole coefficients, which are a priori the most challenging to recover, are shown in Figure 5 after k = 10 iterations in the case of an underlying isotropic flux of cosmic rays. The reconstructed histograms are observed to be well-described by Gaussian functions centered on zero and with a dispersion following indeed Equation (11).

Zoom In Zoom Out Reset image size

With exactly the same ingredients, the simulations can be used to test the procedure with an underlying anisotropic flux of cosmic rays, chosen here such that Φ(n)∝1 + 0.1 Y₁₀(n) + 0.1 Y₂₀(n). Results of the Monte-Carlo simulations are shown in Figure 6 for the specific a₁₀ and a₂₀ coefficients. Again, the reconstructed histograms are observed to be well described by Gaussian functions with parameters following the expectations.

Zoom In Zoom Out Reset image size

Note that in practice, all results presented here are found to be stable as soon as the number of iterations is k = 4.

Formally, the implementation of the cross-calibration procedure is not limited to the choice of the whole overlapping declination band for the integration range ΔΩ in previous equations. The choice of the whole common band turns out to be, however, optimal in terms of resolution in b. A restriction of the common declination band to, for instance, [ − 10°, 10°] would lead to a resolution in b of ≃ 5%; while the use of the whole sky would not bring any improvement for resolving better b. In next sections, the cross-calibration procedure is thus applied to the joint data set by using the whole overlapping region for ΔΩ.

All analyses reported in this section are based on a joint data set consisting of events with energies in excess of roughly 10¹⁹ eV in terms of the energy scale used at the Telescope Array by evaluating in the Auger data set the energy threshold which guarantees equal fluxes for both experiments. We are thus left here with 2130 events (795 in the common band) above 10¹⁹ eV from the Telescope Array and 11,087 (3435 in the common band) above 8.5 × 10¹⁸ eV from the Pierre Auger Observatory. The arrival directions are shown in Figure 7 in equatorial coordinates using the Mollweide projection. Auger data can be seen as the red points in the Southern hemisphere, while Telescope Array ones are shown as the black crosses in the Northern hemisphere.

Zoom In Zoom Out Reset image size

The methodology presented in the previous section allows us to estimate the multipole coefficients and to perform a rich series of anisotropy searches by taking profit of the great advantage offered by the full-sky coverage. After iterations, the coefficient b is b = 1.011. Choosing to use nominal energies to build the joint data set would lead to a different value for b (0.755) due to the different statistics in the Auger data set (8259 events in total instead of 11,087), but it will be shown in next section that this choice impacts the physics results to only a small extent.

The normalization convention of the multipole moments used hereafter is chosen so that the a_ℓm coefficients measure the relative deviation with respect to the whole contribution of the monopole (i.e., the a_ℓm coefficients are redefined such that $a_{\ell m} \rightarrow \sqrt{4\pi }a_{\ell m}/a_{00}$).

5.1. Multipolar Analysis

The dipole, quadrupole, and octupole moments as derived from the iterative procedure are reported in Table 1 in equatorial coordinates together with their associated uncertainties calculated from Equation (11). None of these multipole coefficients stands out as being significantly above the noise level.

Table 1. First Low-order Multipolar Moments and Their Uncertainties (in Equatorial Coordinates)

ℓ	m	aℓm		ℓ	m	aℓm		ℓ	m	aℓm
									−3	−0.022 ± 0.034
					−2	0.038 ± 0.035			−2	0.030 ± 0.039
	−1	−0.102 ± 0.036			−1	0.067 ± 0.040			−1	0.067 ± 0.037
1	0	−0.006 ± 0.074		2	0	0.017 ± 0.042		3	0	−0.027 ± 0.040
	1	−0.001 ± 0.036			1	0.004 ± 0.040			1	0.009 ± 0.037
					2	0.040 ± 0.035			2	−0.004 ± 0.039
									3	−0.011 ± 0.034

Download table as:  ASCIITypeset image

The full set of multipole coefficients provides a comprehensive description of the anisotropy patterns that might be present in the data. A significance table for the coefficients up to ℓ = 20, built simply by dividing each estimated coefficient by its corresponding uncertainty, is reported in the left panel of Figure 8. As it can be seen from the contrast scale, significance values between −1 and 1 dominate the picture. Deviations close to −3 and 3 stand at the expected level for isotropy, as shown in the right panel. Hence, overall, the extraction of the multipole coefficients does not provide any evidence for anisotropy.

Zoom In Zoom Out Reset image size

5.2. Flux and Overdensities/Underdensities Sky Maps

To visualize the result of the multipolar expansion, a flux sky map of the joint data set is displayed using the Mollweide projection in the left panel of Figure 9, in units of km⁻² yr⁻¹ sr⁻¹. This map is drawn in equatorial coordinates. To exhibit structures at intermediate scales, the expansion is truncated here at ℓ = 4. Relative excesses and deficits are clearly visible on a 15% contrast scale.

Zoom In Zoom Out Reset image size

To quantify whether or not some contrasts are statistically compelling, a significance sky map of the overdensities/underdensities obtained in circular windows of radius 15° is shown in the right panel. The choice of the 15° angular scale is well suited to exhibit structures at scales that can be captured by the set of low-order multipoles up to ℓ = 4. Significances are calculated using the widely used Li and Ma estimator (Li & Ma 1983), S, which was designed to account for both the fluctuations of the background and of an eventual signal in any angular region searched,

$\begin{eqnarray} S&=&\pm \sqrt{2}\left[N_{\rm on}\ln {\frac{(1+\alpha _{\rm LM})N_{\rm on}}{\alpha _{\rm LM}(N_{\rm on}+N_{\rm off})}}\right.\nonumber\\ &&\left.+\,N_{\rm off}\ln {\frac{(1+\alpha _{\rm LM})N_{\rm off}}{N_{\rm on}+N_{\rm off}}}\right]^{1/2}, \end{eqnarray} \tag{ 15 }$

with N_on the observed number of events in the angular region searched, and N_off the total number of events. The sign of S is chosen positive in case of excesses and negative in case of deficits. On the other hand, since the background estimation is not based here on any on/off procedure but can be instead determined from the integration of the directional exposure in the angular region searched, the α_LM parameter expressing the expected ratio of the count numbers between the angular region searched and any background region is taken here as

$\begin{equation} \alpha _{\rm LM}(\mathbf {n})=\frac{\int {d}\mathbf {n}^{\prime }\,\bar{\omega }(\mathbf {n}^{\prime })\,f(\mathbf {n},\mathbf {n}^{\prime })}{ \int {d}\mathbf {n}\,\bar{\omega }(\mathbf {n})}, \end{equation} \tag{ 16 }$

with f the top-hat filter function at the angular scale of interest. In the absence of signal, the variable S is expected to be nearly normally distributed. Hence, for positive (negative) values, S (−S) can be interpreted as the number of standard deviations of any excess (deficit) in the sky.

Overall, overdensities and underdensities obtained in circular windows of radius 15° are well reproduced by the multipolar expansion. Contrasts are not identical in all regions of the sky due to the non-uniform coverage (high flux values in low-exposed regions can lead to overdensities less significant than lower flux values in higher exposed regions, and vice versa), but the overall pattern looks similar between the two maps. From the significance contrast scale in the right panel, it is clear that there is no overdensity or underdensity standing above the 3 standard deviation level. The distribution of significances turns out to be compatible with that expected from fluctuations of an isotropic distribution.

5.3. Dipole and Quadrupole Moments

As outlined in the Introduction, although the full set of spherical harmonic moments is needed to characterize any departure from isotropy at any scale, the dipole and quadrupole moments are of special interest. For that reason, special emphasis is given here to these low-order moments, in terms of a more traditional and geometric representation than the raw result of the multipole moments.

The dipole moment can be fully characterized by a vector with an amplitude r and the two angles {δ_d, α_d} of the unit vector d. The quadrupole, on the other hand, can be fully determined by two independent amplitudes {λ₊, λ₋}, two angles $\lbrace \delta _{{\rm q}_+},\alpha _{{\rm q}_+}\rbrace$ defining the orientation of a unit vector q₊, and one additional angle $\alpha _{q_-}$ defining the directions of another unit vector q₋ in the orthogonal plane to q₊. The full description is completed by means of a third unit vector q₀, orthogonal to both q₊ and q₋, and with a corresponding amplitude such that the traceless condition λ₊ + λ₀ + λ₋ = 0 is satisfied. The estimation of the amplitudes and angles of the unit vectors from the estimated spherical harmonic moments is straightforward (see the Appendix). The parameterization of the low-order moments of the flux is then written in a convenient and intuitive way as

$\begin{eqnarray} \Phi (\mathbf {n})&=&\frac{\Phi _0}{4\pi }(1+r\,\mathbf {d}\cdot \mathbf {n}+\lambda _+(\mathbf {q}_+\cdot \mathbf {n})^2\nonumber\\ &&+\lambda _0(\mathbf {q}_0\cdot \mathbf {n})^2+\lambda _-(\mathbf {q}_-\cdot \mathbf {n})^2+\cdots ). \end{eqnarray} \tag{ 17 }$

The values of the estimated amplitudes and angles are given in Table 2 with their associated uncertainties. The distributions of amplitudes obtained from statistical fluctuations of simulated isotropic samples are shown in Figure 10. The superimposed arrows, indicating the measured values, are clearly seen to stand within high probable ranges of amplitudes expected from isotropy.

Zoom In Zoom Out Reset image size

Table 2. Amplitudes and Angles of the Dipole Vector and Quadrupole Tensor

	Amplitude	δ	α	l	b
	(%)	(°)	(°)	(°)	(°)
d	5.0 ± 1.8	3 ± 30	89 ± 22	204	−10
q+	4.0 ± 1.8	−42 ± 41	46 ± 69	260	58
q−	−5.3 ± 2.0	28 ± 22	106 ± 76	154	25
q0	1.3 ± 1.6	34 ± 29	354 ± 72	113	−5

Download table as:  ASCIITypeset image

The dipole parameters, namely, the amplitude, declination and phase, are observed to be compatible with previous reports from both experiments (The Pierre Auger Collaboration 2011a, 2012, 2013; Tinyakov et al. 2012). It is worth noting this for the phase α_d of the dipole vector d: a consistency of phases in adjacent energy intervals was also pointed out in the Auger data (The Pierre Auger Collaboration 2011a, 2012). Given that a consistency of phases is expected to manifest with a smaller number of events than those required for the detection of significant amplitudes, continued scrutiny of future data will provide evidences of whether this phase consistency in both hemispheres is indicative of a real anisotropy or not.

5.4. Power Spectrum

The angular power spectrum C_ℓ is a coordinate-independent quantity, defined as the average |a_ℓm|² as a function of ℓ,

$\begin{equation} C_\ell =\frac{1}{2\ell +1}\sum _{m=-\ell }^\ell |a_{\ell m}|^2. \end{equation} \tag{ 18 }$

In the same way as the multipole coefficients, any significant anisotropy of the angular distribution over scales near 1/ℓ radians would be captured in a non-zero power in the mode ℓ. Although the exhaustive information of the distribution of arrival directions is encoded in the full set of multipole coefficients, the characterization of any important overall property of the anisotropy is hard to handle in a summary plot from this set of coefficients. Conversely, the angular power spectrum does provide such a summary plot. In addition, it is possible that for some fixed mode numbers ℓ, all individual a_ℓm coefficients do not stand above the background noise but meanwhile do so once summed quadratically.

From the set of estimated coefficients $\bar{a}_{\ell m}$, the measured power spectrum is shown in Figure 11. The gray band stands for the rms of power around the mean values expected from an isotropic distribution, while the solid line stands for the 99% confidence level upper bounds that would result from fluctuations of an isotropic distribution. Overall, no significant deviation from isotropy is found from this study.

Zoom In Zoom Out Reset image size

There are uncertainties in choosing the energy scales to be used when building the joint data set, and/or in correcting the energy estimator for instrumental effects, and these propagate into systematic uncertainties in the measured anisotropy parameters. In this section, we choose to use the angular power spectrum as a relevant proxy to probe the size of the systematic effects investigated below.

Even though the cross-calibration of energies can be considered as a reasonable starting point for building the joint data set, it is not a necessary input for the iterative procedure described in Section 4. Using instead the nominal energies of each experiment only results in a different balance between the two nominal exposures obtained by requiring equal fluxes in the common band through Equation (12). Then, the final anisotropy parameters necessarily differ, but only slightly. This is evidenced in Figure 12, where the points (and their statistical uncertainties) of the power spectrum obtained previously are shown as the gray bands for each moment ℓ. The power spectrum obtained using nominal energies is shown as the filled circles. The picture is in global agreement, with only a few points standing at most one standard deviation away.

Zoom In Zoom Out Reset image size

As already mentioned in Section 2.3, some distortions are imprinted in the event rate as a function of local angles by the influence of the geomagnetic field on the development of the showers. The importance of this effect depends on the weight of the muonic component to the signal size. At the Pierre Auger Observatory, this effect has been shown to induce significant variations of the event rate in declination as soon as the total number of events analyzed is of the order of 30,000 (The Pierre Auger Collaboration 2011b) if the corrections of the energy estimator discussed in Section 2.3 are not applied. Given the current statistics above 10¹⁹ eV, the distortions are however expected to be marginal for the specific analysis reported here. This is evidenced by the power spectrum shown as the open squares in Figure 12 obtained without applying the corrections to the (cross-calibrated) Auger data set. All points are indeed observed to be within the statistical uncertainties of the estimate shown in Figure 11. Given the smaller statistics available in the data set of the Telescope Array, and given that the size of this geomagnetic effect is expected to be smaller due to the lower weight of the muonic component to the signal size with scintillators compared to water-Cherenkov detectors, this provides support that the absence of energy corrections in the data set of the Telescope Array does not impact on the results presented in this report.

Note that the spread of the 99% confidence level line in Figure 12 stands for the slightly different statistics which result when using nominal or cross-calibrated energies to select all events above 10¹⁹ eV.

In this work, an entire mapping of the celestial sphere has been presented by combining data sets recorded at the Telescope Array and the Pierre Auger Observatory above 10¹⁹ eV. The unavoidable systematic uncertainty in the relative exposures has been treated by designing a cross-calibration procedure relying on the common region of the sky covered by both experiments. This cross-calibration procedure makes it possible to use the powerful multipolar analysis method for characterizing the sky map of ultra-high energy cosmic rays. Throughout the series of anisotropy searches performed, no significant deviation from isotropy could be captured at any angular scale.

From the multipolar coefficient measurements performed in this work, upper limits can be derived for any kind of pattern. The ones obtained at 99% confidence level on the dipole and quadrupole amplitudes are shown in Figure 13 as a function of the direction in the sky, in Equatorial coordinates. These upper limits have been obtained by searching for the smallest values of dipole amplitude oriented in each direction d and quadrupole amplitude oriented in each direction q₊ guaranteeing that the reconstructed amplitudes in simulated data sets are larger or equal to the ones obtained for real data in 99% of the simulations. The different sensitivities for each direction are caused by the different resolutions for each reconstructed multipolar coefficient. Note that the upper limits on the quadrupole amplitude pertain to a symmetric quadrupole only (that is, a quadrupole with amplitudes such that λ₋ = λ₀ = −λ₊/2) to keep the number of studied variables manageable.

Zoom In Zoom Out Reset image size

For the first time, the upper limits on the dipole amplitude reported in Figure 13 do not rely on any assumption on the underlying flux of cosmic rays thanks to the full-sky coverage achieved in this joint study. With partial-sky coverage, similar sensitivity could be obtained in this energy range but assuming a pure dipolar flux (The Pierre Auger Collaboration 2012). In addition, the sensitivity on the quadrupole amplitude (and to higher order multipoles as well) is the best ever obtained thanks, also, to the full-sky coverage.

The cross-calibration procedure designed in this study pertains to any combined data sets from different observatories showing an overlap in their respective directional exposure functions and covering the whole sky once combined. It is conceivable to apply it in an energy range where the detection efficiency is not saturated. Then, future joint studies between the Telescope Array and the Pierre Auger collaborations will allow further characterization of the arrival direction distributions down to ≃ 10¹⁸ eV.

The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe. We are very grateful to the following agencies and organizations for financial support: Comisión Nacional de Energía Atómica, Fundación Antorchas, Gobierno De La Provincia de Mendoza, Municipalidad de Malargüe, NDM Holdings, and Valle Las Leñas, in gratitude for their continuing cooperation over land access, Argentina; the Australian Research Council; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (FINEP), Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ), São Paulo Research Foundation (FAPESP) grants # 2010/07359-6 and # 1999/05404-3, Ministério de Ciência e Tecnologia (MCT), Brazil; MSMT-CR LG13007, 7AMB14AR005, CZ.1.05/2.1.00/03.0058, and the Czech Science Foundation grant 14-17501S, Czech Republic; Centre de Calcul IN2P3/CNRS, Centre National de la Recherche Scientifique (CNRS), Conseil Régional Ile-de-France, Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS), Département Sciences de l'Univers (SDU-INSU/CNRS), Institut Lagrange de Paris, ILP LABEX ANR-10-LABX-63, within the Investissements d'Avenir Programme ANR-11-IDEX-0004-02, France; Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Finanzministerium Baden-Württemberg, Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF), Ministerium für Wissenschaft und Forschung, Nordrhein Westfalen, Ministerium für Wissenschaft, Forschung und Kunst, Baden-Württemberg, Germany; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR), Gran Sasso Center for Astroparticle Physics (CFA), CETEMPS Center of Excellence, Italy; Consejo Nacional de Ciencia y Tecnología (CONACYT), Mexico; Ministerie van Onderwijs, Cultuur en Wetenschap, Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Stichting voor Fundamenteel Onderzoek der Materie (FOM), The Netherlands; National Centre for Research and Development, grant Nos. ERA-NET-ASPERA/01/11 and ERA-NET-ASPERA/02/11, National Science Centre, grant Nos. 2013/08/M/ST9/00322, 2013/08/M/ST9/00728, and HARMONIA 5 - 2013/10/M/ST9/00062, Poland; Portuguese national funds and FEDER funds within COMPETE - Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia, Portugal; Romanian Authority for Scientific Research ANCS, CNDI-UEFISCDI partnership projects nr.20/2012 and nr.194/2012, project nr.1/ASPERA2/2012 ERA-NET, PN-II-RU-PD-2011-3-0145-17, and PN-II-RU-PD-2011-3-0062, the Minister of National Education, Programme for research—Space Technology and Advanced Research—STAR, project number 83/2013, Romania; Slovenian Research Agency, Slovenia; Comunidad de Madrid, FEDER funds, Ministerio de Educación y Ciencia, Xunta de Galicia, European Community 7th Framework Program, grant No. FP7-PEOPLE-2012-IEF-328826, Spain; Science and Technology Facilities Council, United Kingdom; Department of Energy, contract No. DE-AC02-07CH11359, DE-FR02-04ER41300, and DE-FG02-99ER41107, National Science Foundation, grant No. 0450696, The Grainger Foundation, USA; NAFOSTED, Vietnam; Marie Curie-IRSES/EPLANET, European Particle Physics Latin American Network, European Union 7th Framework Program, grant No. PIRSES-2009-GA-246806; and UNESCO.

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science through Grants-in-Aids for Scientific Research on Specially Promoted Research (21000002) "Extreme Phenomena in the Universe Explored by Highest Energy Cosmic Rays" and for Scientific Research (19104006), and the Inter-University Research Program of the Institute for Cosmic Ray Research; by the U.S. National Science Foundation awards PHY-0307098, PHY-0601915, PHY-0649681, PHY-0703893, PHY-0758342, PHY-0848320, PHY-1069280, and PHY-1069286; by the National Research Foundation of Korea (2007-0093860, R32-10130, 2012R1A1A2008381, 2013004883); by the Russian Academy of Sciences, RFBR grants 11-02-01528a and 13-02-01311a (INR), IISN project No. 4.4509.10 and Belgian Science Policy under IUAP VII/37 (ULB). The foundations of Dr. Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles and the George S. and Dolores Dore Eccles all helped with generous donations. The State of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management, and the U.S. Air Force. We also wish to thank the people and the officials of Millard County, Utah for their steadfast and warm support. We gratefully acknowledge the contributions from the technical staffs of our home institutions. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged.

We provide in this appendix the transformation rules between the multipole coefficients and the parameters of the dipole vector and the quadrupole tensor. The multipole coefficients are assumed to be calculated from arrival directions expressed in equatorial coordinates. The Cartesian components of the dipole vector d are related to the a_1m coefficients through

$\begin{equation} d_x=\frac{\sqrt{3}}{a_{00}}\,a_{11},\quad d_y=\frac{\sqrt{3}}{a_{00}}\,a_{1-1},\quad d_z=\frac{\sqrt{3}}{a_{00}}\,a_{10}. \end{equation} \tag{ A1 }$

The amplitude d and directions δ_d and α_d are then obtained by

$\begin{equation} d=\sqrt{d_x^2+d_y^2+d_z^2},\quad \delta _d=\arcsin {d_z},\quad \alpha _d=\arctan {d_y/d_x}. \end{equation} \tag{ A2 }$

The quadrupole can be described by a second order tensor Q such that the flux can be expressed as

$\begin{equation} \Phi (\mathbf {n})=\frac{\Phi _0}{4\pi }\left(1+r\,\mathbf {d}\cdot \mathbf {n}+\frac{1}{2} {\sum _{i,j}Q_{ij}\,n_i\,n_j}\right). \end{equation} \tag{ A3 }$

Q is a traceless and symmetric tensor. Its five independent components are determined from the a_2m by

$\begin{eqnarray} Q_{xx} &= \frac{\sqrt{5}}{a_{00}}\,(\sqrt{3}a_{22}-a_{20}), \end{eqnarray} \tag{ A4 }$

$\begin{eqnarray} Q_{xy} &= \frac{\sqrt{15}}{a_{00}}\,a_{2-2}, \end{eqnarray} \tag{ A5 }$

$\begin{eqnarray} Q_{xz} &= -\frac{\sqrt{15}}{a_{00}}\,a_{21}, \end{eqnarray} \tag{ A6 }$

$\begin{eqnarray} Q_{yy} &= \frac{\sqrt{5}}{a_{00}}\,(\sqrt{3}a_{22}+a_{20}), \end{eqnarray} \tag{ A7 }$

$\begin{eqnarray} Q_{yz} &= -\frac{\sqrt{15}}{a_{00}}\,a_{2-1}. \end{eqnarray} \tag{ A8 }$

The other components are obtained by symmetry and from the traceless condition (that is, Q_zz = −Q_xx − Q_yy). The amplitudes λ_{±, 0} are then obtained as the eigenvalues of Q and the vectors q_{±, 0} as the corresponding eigenvectors.