Beyond the Galaxy: UHECR results from the Pierre Auger Observatory and the Telescope Array

The beginning of the 21st century has brought much progress into the field of Ultra High Energy Cosmic Rays particularly through the completion of the Pierre Auger Observatory first and the Telescope Array later on. The success of these detectors follows from their hybrid character combining the advantages of the fluorescence technique for energy determination with the larger exposure provided by an array of particle detectors. In this article we review the important contributions made by both experiments concerning the spectrum, anisotropy searches and composition studies. We discuss how these measurements have contributed to progress in the field, partially solving long standing puzzles, and how they have also led to new challenges that are likely to constitute the driving force of the field in the immediate future.


Introduction
The existence of ultra high energy cosmic rays reaching energies in the vicinity of 10 20 eV (100 EeV) has been known for about half a century now [1]. These extremely energetic particles are uniquely interesting for multiple reasons both in astrophysics and in particle physics. It is difficult to envisage astrophysical objects capable of accelerating particles to such extreme energies or alternative mechanisms of production. Because of their large rigidity, they approximately keep the direction with which they are produced over paths greatly exceeding galactic dimensions and, if they are of extragalactic origin, which seems to be most plausible, they must have short ranges on cosmological scales due to interactions with the Cosmic Microwave Background (CMB) and with the intergalactic magnetic fields. Their collisions with the CMB must produce fluxes of neutrinos in the EeV range which could be detected, while their collisions with atmospheric nuclei have center of mass energies that exceed those achieved in particle accelerators by several orders of magnitude. Not surprisingly they have attracted much attention and many efforts to observe them with increasing statistics and precision since the early days. In spite of enormous progress there are still crucial issues which are unresolved. Finding answers to the nature and origin of these particles has a high priority since they are likely to forge the future of the field.
The two largest and most precise detectors to date are the Pierre Auger Observatory in Argentina (Mendoza) and the Telescope Array in the US (Utah). Both detectors are hybrid, combining an array of particle detectors to sample extensive air showers when they reach the ground and telescopes to collect the fluorescence light of the atmospheric nitrogen excited by billions of shower particles crossing it. This combination of techniques optimizes their performance and keeps the need to use simulations to a bare minimum.
The two experimental facilities have accumulated large exposures during years of operation and their data have provided measurements of the spectrum, the arrival directions and the shower development with unprecedented precision, leading to significant advances in the field such as the establishment of the suppression of the flux at the highest energies. The interpretation of these measurements is however far from settled, particularly concerning the arrival directions and the primary composition since both experiments arrive at different conclusions. I here review the progress made in the measurement of the spectrum, the study of anisotropies, and the primary mass composition of these enigmatic particles.
2. The success of the hybrid approach 2.1. Auger Observatory and Telescope Array Both the Pierre Auger Observatory and the Telescope Array are based on similar detector concepts. They are both located at roughly similar latitudes respectively in the southern and northern hemispheres and at an average altitude of 1,400 m. The main characteristics of these detectors are summarized in table 1 and further details can be found elsewhere [2,3,4]. The most fundamental difference lies in their surface detectors (SD) which are based on different detector concepts. The particle detectors in the Auger SD are cylindrical tanks of 10 m 2 surface and 1.2 m height filled with purified water, with three photomultiplier tubes (PMT) to detect the Cherenkov light of particles in the shower front. In the TA they consist of two 3 m 2 slabs of plastic scintillator on top of one another which give light pulses also read by PMTs. The former are sensitive to the track of the charged particles and the latter count charged particles as they cross. The water tanks are relatively much more sensitive to shower muons which usually traverse the tank from wall to wall while the counts in the TA detectors are dominated by electrons and positrons that typically outnumber the muons in the shower front. The detectors of the Pierre Auger Observatory are well suited to detect highly inclined showers, increasing the exposure and the sky coverage. The analysis of inclined showers has also been used for neutrino searches and to establish the muon content of the showers.
The Auger Observatory extends over an area of 3,000 km 2 while the TA covers about 700 km 2 . The fluorescence telescopes are located on the perimeter of the two observatories to monitor practically the whole atmospheric volume just above the surface arrays. They have photomultiplier cameras in the focal plane which measure the arrival time, direction and intensity of the emitted light as the shower develops. Some details like the telescope number, size, pixeling and viewing angles can be compared in the aforementioned table. The Pierre Auger Observatory contains a smaller area of 23.5 km 2 with stations separated by 750 m (750 m array) which can be combined with three additional telescopes pointing at higher elevations (HEAT) for lower energy measurements. Similarly the Telescope Array has two sub-arrays of 46 and 35 stations separated by about 600 m and 400 m over a surface of ∼20 km 2 , together with ten telescopes covering from 31 • to 59 • of elevation (TALE). These detectors are designed to measure the spectrum for energies above 30 PeV [5]. The Pierre Auger Observatory is run by an international collaboration involving over 500 scientists from more than 100 institutions in 19 countries. The Telescope Array is run by about 150 scientists from 27 institutions in 5 countries.
The combination of two detection techniques which are well established has multiple advantages. The fluorescence detector (FD) allows a calorimetric measurement of the energy deposition in the atmosphere providing a reliable determination of the shower energy. In addition it also allows a measurement of the depth of shower maximum which gives a good handle on primary composition. On the other hand the surface detector (SD) has a very well defined exposure, related to the geometrical area of the surface array, and a duty cycle of nearly 100%, crucial to improve the statistical precision. Furthermore the redundancy introduced by the two  [7], amounts to first establishing the arrival direction of the shower using timing information and then converting the measured light signal to the energy deposition in the atmosphere after accounting for absorption and dispersion. Sophisticated monitoring systems of the atmosphere are implemented in both observatories to minimize the uncertainties in the calculation of the emitted light. If the shower is viewed from a single site, the uncertainty in the arrival direction can be quite large. Both collaborations use the hybrid potential to get around this problem, combining timing information of at least one SD station with that of the signals in the pixels of the camera. The amount of emitted light is assumed to be proportional to energy deposition and the relation between particle numbers and emitted light is given by the fluorescence yield which has been established experimentally. Unfortunately both collaborations use different parameterizations of the fluorescence yield, which complicate comparisons between results.

2.2.2.
Surface array for zenith angles below 60 • The reconstruction of events with the SD is very similar in the two observatories. The first step consists on establishing the arrival direction using the arrival times of the signals in the particle detectors. Once the direction is known, the pattern of signals on the ground is fitted to obtain the position of the shower core and the shower size. For showers with zenith angles, θ, below 60 • the signal pattern can be assumed to have approximate circular symmetry around the shower axis so that it is parameterized as a function of a single variable, r, the distance to shower axis. For a fixed mass composition and arrival direction the shower size scales with energy, and it is characterized by the value of the fitted signal, S, at a reference distance, chosen to minimize fluctuations. These parameters are S(1000) at 1000 m for the Auger data with θ < 60 • , S(800) at 800 m for the TA data, restricted to θ < 45 • , and S(450) at 450 m for the data of the Auger 750 m array with θ < 55 • . The TA Collaboration compares measured values of S(800) to simulations to obtain a first estimate of shower energy. The conversion is made with a table obtained from detailed simulations of extensive air showers and the signals they produce on the SD. These simulations have been made assuming constant proton composition, an energy spectrum compatible to earlier measurements, and using the QGSJETII-03 [8] hadronic interaction model. The Auger Collaboration however first finds a quantity that scales with energy and which is independent of zenith angle to use it as an energy estimator. For the regular (750 m) array this quantity is chosen to be S 38 (S 35 ), the value of S(1000) (S(450)) that the shower would have had if it had the same energy but arrived at the median zenith angle of 38 • (35 • ). Assuming isotropy of the arrival directions the conversion is made with the constant intensity cut method, avoiding the use of simulation [9].
2.2.3. Surface array for zenith angles exceeding 60 • For high zenith angles the geomagnetic field breaks the cylindrical symmetry and a specific method has been devised to fit the expected patterns on the ground to the data [10]. For a given arrival direction the number density of muons has a two-dimensional pattern [11] approximately universal in shape, independently of shower energy, primary composition and hadronic model used [12].
For the Auger data the distributions used for reference correspond to the average values of 10 EeV-proton showers simulated with QGSJETII-03 which have been conveniently parameterized [13]. The fitted normalization, N 19 , is a size parameter which is by construction independent of arrival direction and thus can be used as an energy estimator.  Correlation between the FD energy and SD energy for data of the Telescope Array which can be simultaneously reconstructed with the FD and SD techniques. The SD energy, obtained by directly comparing S(800) to simulations, has been rescaled by a constant factor of 1.27.

2.2.4.
Energy calibration The showers detected with the surface detector of the two observatories have been calibrated with their fluorescence detectors. Both collaborations take advantage of the hybrid measurements using a sub-sample of events well reconstructed with both the SD and FD techniques. This allows them to translate the more precise measurement of the energy obtained with the FD, which is however limited in duty cycle to moonless clear nights, to the larger statistical sample obtained with the SD array. The Pierre Auger Collaboration calculates the correlation between the energy estimator and the FD energy [14]. A relation between FD energy and the value of each of the three estimators is obtained as illustrated in Fig.1. The TA has found that rescaling the derived SD energy with a constant factor of 1.27 ensures a good matching to the measured FD energy as illustrated in Fig. 2. The plot is rather similar to the calibration plots of the Pierre Auger Collaboration but we note that the TA procedure does not allow a deviation from linear scaling between the measured energy and the inferred energy from simulations. Such deviations should however be expected in the case that there was a significant change of composition in the explored energy range. Shower simulations show that when composition is changed for a given energy and arrival direction, the average position of shower maximum the average number of muons changes. Consequently the average signal at a given distance to shower axis differs depending on mass composition.

The shape of the UHECR energy spectrum
The spectrum is calculated selecting data that fall inside a given region of the detector, and dividing the energy distribution by the corresponding exposure. If the showers have sufficient energy, the trigger efficiency of the surface detector becomes close to 100% and the exposure is independent of energy and simply obtained by integrating over the surface area, time and solid angle. The SD of the Pierre Auger Observatory becomes fully efficient above 3 EeV for events with zenith angle, θ, below 60 • . For inclined events with θ exceeding 60 • the energy above which the efficiency is approximately 100% is 4 EeV. For events detected with the 750 m array full efficiency is reached at energies greater than 0.3 EeV. The TA Collaboration includes the efficiency in the calculation of the exposure and measures the spectrum for energies above 1.6 EeV. At the lower energies the detection efficiency is of order 10%, and uncertainties in the efficiency directly translate into the flux determination. Full efficiency is achieved at the TA for energies above about 10 EeV. Three independent fluxes have been obtained with the Auger Observatory at the highest energies [15,16]. The data collected up to 31st December 2012 with θ < 60 • give the largest exposure corresponding to 82,318 events of energies above 3 EeV. The highly inclined data (62 • < θ < 80 • ) with about a third of the exposure correspond to 11,074 events of energy above 4 EeV and the hybrid data add up to 11,155 events above 1 EeV during the same period. The results of these three fluxes and the flux obtained at lower energies with the 750 m array are compared in Fig. 3. The Telescope Array with about a tenth of the total Auger SD exposure has collected 14,747 events above 1.6 EeV during the period between 11 May 2008 and 20 May 2012.
The results of the spectrum obtained with the TA data are compared in Fig. 4 to those obtained by the Pierre Auger Collaboration after combining the four spectra obtained with hybrid data and with each of the three sets of SD data (θ < 60 • , very inclined and 750 m array). They provide unambiguous evidence of flux suppression at the highest energies. In addition a break in the spectrum is clearly appreciated, the ankle, at about 5 EeV. Above this energy, the flux is harder. The agreement between the spectra of the two observatories for energies below 40 EeV is excellent, well within the systematic uncertainties. The ankle energy is found at 5.25 ± 0.27 EeV and 4.6 ± 0.3 EeV respectively for the Auger and TA data and the spectral indices at either side are fully compatible. The difference between the fluxes is of the same order as that obtained when the convention of the fluorescence yield from one experiment is changed to that used by the other.
At the highest energies, the differences between the two fluxes are much more significant leading to different interpretations. The energy at which the integrated flux is half of that  Figure 3. Comparison of the four spectra obtained with the Auger Observatory. The spectrum obtained with the SD for events with θ < 60 • is shown with black squares the inclined spectrum (62 • < θ < 80 • ) is plotted with red upright triangles and that obtained in the 750 array is indicated with blue circles. The spectrum obtained with hybrid data is indicated with green triangles. Only statistical uncertainties are indicated with the error bars.
expected by extrapolating the spectral index found above the ankle is usually denoted as E 1/2 . The Pierre Auger Collaboration finds E 1/2 = 43 ± 1 EeV to be compared with 60 EeV for the Telescope Array. Besides the fluorescence yield, we note differences in the treatment of the missing energy correction (energy that is not deposited in the atmosphere, mainly attributed to neutrinos and muons) in the procedure used for calibrating the energy to the fluorescence measurement and note that the two experiments are looking at opossite parts of the sky. It is possible that composition changes and/or anisotropies at these energies contribute to these differences. The shape of the spectrum alone can give us information about possible sources and locations of UHECR. It is tempting to favor a pure proton scenario in which both the suppression and the ankle are features naturally resulting from interactions with the CMB, as the TA Collaboration claims. Models in which protons are produced by isotropically distributed sources emitting protons with free parameters such as the maximum accelerating energy, the spectral index of the emission at the sources and the cosmological evolution of the sources have been successfully fit to the spectral data [17]. The same exercise has been done by the Pierre Auger Collaboration which has also found good fits to the data. However these fits reveal that the maximum accelerating  energy of these sources must be about 100 EeV. This implies that the observed suppression at the highest energies is an intrinsic property of the sources and no longer the effect of the interactions with the CMB alone. A pure proton composition is not consistent with other data from the Auger Observatory as will be discussed below. Other scenarios with pure iron or mixed composition have also been fitted to the spectrum measured by the Pierre Auger Collaboration. In these scenarios the ankle is interpreted as the transition from galactic to extragalactic origin [18] and ceases to be a natural consequence of particle interactions unless the spectrum is dominated by protons at the highest energies.

Hints of anisotropies
The arrival directions provide an important handle on fundamental issues relating to the origin and nature of UHECR. Magnetic deviation can be sufficiently small if these particles have large rigidities and they are not produced too far away. This would be the case if the transition from galactic to extragalactic origin is responsible for the ankle feature. Large-scale anisotropies would be expected in the EeV range, signalling particles that escape from the galaxy. The expected patterns are however not easily predictable due to the uncertainty in magnetic fields and the unknown distribution of sources. Naturally, the patterns would be much stronger in case the composition is light, in agreement with most EeV data. Even in the case that these cosmic rays are extragalactic and are highly isotropic, our bulk motions relative to the rest mass frame of the cosmic rays would be expected to provide some large-scale anisotropy signal at the level of about 0.6% due to our relative motions with respect to the rest frame of the CMB radiation [19]. Above the 50 EeV region, interactions with the background radiation fields restrict the range of proton and iron particles to distances of order 200 Mpc or less and this becomes even lower for particles of intermediate charge [20]. Intergalactic magnetic fields are often assumed to be coherent over regions of ∼1 Mpc with intensity below ∼1 nG which leads to a few degrees deflection for protons over these distances [21]. As the matter distribution is not homogeneous in our vicinity at scales of order 100 Mpc it is again reasonable to expect anisotropies, particularly if the particles are protons. The attractive possibility of doing astronomy with charged particles requires light mass composition at the highest energies. Composition measurements are needed. Both the Pierre Auger and the TA Collaborations have extensive programs to search for and measure anisotropies of UHECR. These include auto-correlation functions, correlation with existing catalogs, searches for excess regions in the sky, for point sources, for large scale anisotropies, as well as correlations with other experiments.

Correlations
The first anisotropy signal was claimed by the Pierre Auger Collaboration from the analysis of correlations of arrival directions of UHECR with catalogues of nearby objects. A correlation was reported for events with energies above 57 EeV with objects of redshift below z < 0.018 in the Véron-Cetty and Véron catalogue of Active Galactic Nuclei (AGN), assuming a maximum of 3.1 • directionality loss [22]. The signal was clear and strong at first but weakened to a 38 +7 −6 % correlation when the data increased by a factor of 3 [23] and more recently to 28.1 +3.8 −3.6 % [24] with over ten years of data, to be compared 21% expected from isotropy. The TA Collaboration has also searched for such correlations and found 11 events correlating when only 5.9 are expected from isotropy, what represents a chance probability of 2% [25]. The results of this search were recently reported to be also consistent with the signal reported by the Pierre Auger Collaboration, with a combined chance probability to be due to an isotropic distribution of p = 10 −3 [26]. Cross correlations with other flux limited catalogs of nearby extragalactic objects have also been scanned. The corresponding probabilities for the observed correlations to arise from isotropy are at the percent level. These probabilities are actually heavily penalized for the multiple scans in UHECR energy threshold, maximum object distance, random angular spreading of directions and catalogues.
Analysis for directions with excess number of events over background has been used to search for anisotropies. The search directions can be either arbitrary or preselected according to plausible scenarios or to claims of signals from other experiments. Neither the TA or the Pierre Auger Collaboration have found any significant excess coming from the galactic center [27,28] as claimed by AGASA [29]. The Pierre Auger Collaboration has reported an excess of events in a broad circular region within 18 • of a direction in the vicinity of Centaurus A for events above 57 EeV after a scan. Expectations from isotropy would lead to 3.2 events in this region while 13 events were found [23]. The latest update gives 14 events above 58 EeV in a region within 15 • around the same region [24]. The chance probability of obtaining such an excess from an isotropic distribution has been calculated to be 1.4%, not low enough to claim a signal. It has been noted that this contributes significantly to autocorrelation measurements and to the correlation with different populations of nearby extragalactic objects. More recently, the TA has reported a significant excess of events above the same energy and within a 20 • circle around a direction centered at R.A.(α)=146.7 • , Dec.(δ)=43.2 • in the northern hemisphere [30]. The chance probability first reported for five years of data was equivalent to a 3.4σ effect (posttrial). Including the sixth year increased the chance probability to 4.1σ [30].

Large scale anisotropies
The search for large-scale anisotropies has recently attracted a lot of attention. In the multipolar approach the flux is expanded in spherical harmonics, Y lm (n) to obtain the multipole coefficients of the flux. For a direction-dependent exposure the data are weighted with the inverse of the corresponding relative exposure [31]. Clearly this procedure requires full sky coverage. Otherwise the multipole expansion can be found for the product of the flux times the exposure as measured by the detector. A relation can be found between this expansion and the standard expansion of the flux in spherical harmonics. If the multipole series of the flux alone is assumed to be bounded, the two sets of coefficients are related by a finite matrix multiplication. To obtain the flux coefficients the matrix must be inverted. Unfortunately the uncertainties involved become large when the number of multipoles is high [32].
Alternatively the rotation of the Earth allows one to search for anisotropies because the exposure is approximately invariant in right ascension. The flux is Fourier-expanded in this coordinate in Rayleigh analysis. The first terms can be also obtained taking the difference of fluxes in the East and West sectors. Since both sectors have identical instantaneous exposures, the difference eliminates to first order detector instabilities and weather effects, at the price of a reduced sensitivity. The Fourier coefficients in right ascension can be related to the standard multipole coefficients. When the odd multipoles in the flux expansion are assumed to be reduced to a single dipole (i.e. to l = 1) the first two Fourier terms are exactly a measurement of the dipole projected onto the equatorial plane [33]. In general, this relation holds if the dipole and quadrupole dominate the flux, but it is only approximate.  A modified Rayleigh analysis has been applied to Auger showers of energy exceeding 1 EeV and zenith angles below 60 • including corrections to compensate for small modulations from nonuniform exposure and weather effects. For lower energies, the differential East-West method has been chosen both on the regular and the 750 m arrays [34]. It has been remarked that the amplitudes of the first harmonic in the modulation of the right ascension at a few energy bins lie just beyond the region of 99% C.L. expected from isotropy [35]. Moreover, the dipole phases obtained in adjacent energy bins are consistent with a smooth transition from a dipole direction pointing close to the galactic center at low energies, to one which is rotated by about 180 • at energies above a few EeV. It is expected that in case that some anisotropy exists the first signs would show up precisely in the measured phases that cease to be randomly distributed as expected from statistical noise. A prescription has been established to test these possible hints of anisotropy [34]. The results obtained of the dipole amplitude compared to expectations from isotropy and the dipole phases presented when the prescription was midterm [35], are shown in Figs. 5,6.  Multipole coefficients have also been obtained for energies above 1 EeV, truncating the expansion to l ≤ 2, i.e. assuming that all the anisotropy is reduced to a single dipole or a dipole plus a quadrupole [37]. The 99% C.L. bounds obtained on the dipole and quadrupole amplitudes are at the 2-12 % level. They are illustrated in Figs. 7,8 compared to expectations from generic models in which EeV particles are proton or iron, produced in stationary galactic sources distributed in the galactic disk and emitting isotropically [38]. The regular component of the magnetic field is assumed to have Bisymmetric Spiral Structure with an antisymmetric halo with respect to the galactic plane [39]. Pure iron sources would be marginally allowed and pure protons excluded in this scenario. This is quite restrictive because it will be later shown that iron composition is not at all favored at these energies by either observatory. As a result galactic sources of UHECR are unlikely according to the data. It is not feasible to go further in l in the expansion because the uncertainty becomes too high [37]. This is the consequence of not having full sky coverage. In such case assumptions about the flux itself must be made to obtain the coefficients.

Recent developments
In the time period between the conference and the writing of this article a number of important developments have followed. The combination of the data from the two experiments provides full sky coverage facilitating the direct extraction of the harmonic coefficients for the first time [31]. The energy scales of the two experiments were compared using the declination band in which both detectors overlap as shown in Fig. 9. The results, compatible with isotropy, are limited due   to systematic uncertainties, particularly in the cross calibration of the energy scales. Much more recently, the highly inclined events with zenith angles between 60 • and 80 • and energies above 4 EeV detected with the Pierre Auger Observatory have been for the first time combined with those below 60 • to update large-scale anisotropy searches [36]. Very inclined showers not only represent a further increase of statistics by nearly 30% but they also enhance the sky coverage to 85%, significantly contributing to reduced uncertainties. The exposures of the two sets of data are compared in Fig. 10. Results were obtained for first and second harmonic modulations in both right ascension and in azimuth as well as for multipole coefficients in the hypotheses of either a single dipole or a dipole plus a quadrupole. While the results are compatible with isotropy for energies between 4 and 8 EeV, the raw significance assuming a single dipole is already at the 4σ level for events with E > 8 EeV, thus strengthening the hint of anisotropy. The multipole coefficients assuming a single dipole give a magnitude of 0.073 ± 0.015 in the direction (δ, α) = (95 • ± 13 • , −39 • ± 13 • ). Assuming a dipole plus a quadrupole gives results consistent with those reported above, but the quadrupole coefficients are not significant. Both the dipole strength and the direction are consistent with previous results. A map of the measured flux above 8 EeV smoothed in 45 • windows is shown in Fig 11 resulting in a contrast (maximum relative difference between fluxes) of 21%.

Unravelling UHECR composition
Understanding composition is a very hard task since showers are stochastic processes. Showers of a given energy and direction tend to be rather similar regardless of composition since the shower development is a process by which the initial particle energy is converted to matter, mostly pions and kaons, in successive particle interactions. In addition, all our knowledge about the subtle differences between showers induced by different primaries must be deduced from extensive simulation programs that must extrapolate the properties of particle interactions to center-ofmass energies far exceeding those achieved in accelerator experiments. A number of simulation packets for air showers exist and they can use different models to extrapolate properties of particle interactions such as cross sections, inelasticity and multiplicity. Unfortunately the differences between these extrapolations are subtle and change the shower evolution at levels which are similar to those due to changing the primary composition. As a result, not only has it been notably difficult to extract properties about composition from air shower measurements, but the coupling between an unknown composition and uncertainties in the hadronic models severely limits the conclusions that can be obtained.

Measurements of the average depth of maximum
Composition can be addressed using the depth of the position of the maximum in the number of shower particles, abbreviated as shower maximum, X max , which is measured by the FD. In the simple superposition model, air showers of energy E produced by elements of mass A can be approximated by the superposition of A proton showers of energy E/A. As the shower maximum scales with the logarithm of the traversed depth, showers produced by heavy nuclei are shallower.
In addition the superposition of A showers of lower energy has an X max distribution which is much narrower. In this picture the resulting X max is some kind of average of the maxima of the A single proton showers that make it up, which naturally fluctuates less. Both experiments have made precise measurements of the depth of shower maximum with their FDs that allow a partial reconstruction of shower development. The TA (Auger) Collaboration has reported X max data for energies above 1.6 (0.63) EeV. While the results look at first rather similar, the interpretation made by each collaboration is qualitatively quite different. The TA has interpreted these measurements as compatible with a pure proton composition [40] while the Pierre Auger Collaboration claims that they indicate a composition that changes from being light at 1 EeV to becoming heavier [41].
Part of this discrepancy may lie in the way that the two collaborations address the bias that arises because of the limited field of view of the telescopes and the quality cuts that are required to accurately reconstruct X max . Depending on the zenith angle and impact point of the shower, an FD views a different range of depths as illustrated in Fig. 12 for three shower geometries. A given shower will be selected if the maximum is well within the field of view (shaded region). In general, the distributions will be artificially cut by the field of view and the quality cuts, potentially biasing the result. If there is a mixture of primaries there will also be a composition bias in the selected data.  The Telescope Array Collaboration measures an effective distribution of shower maximum and compares it to simulations which undergo identical selection as the data. Because of the selection bias, the measured average and variance of X max will not match those of the true distribution. The comparison between data and simulation is valid provided the models used accurately describe the ranges of X max . The Pierre Auger Collaboration attempts an unbiased measurement introducing fiducial cuts which reduce the statistics [42]. For each shower that passes the quality cuts the range of values of X max that can be selected must contain that of the bulk of the X max distribution for the corresponding energy characterized by its extremes. These extremes are chosen ensuring that, when they are changed, < X max > remains sufficiently stable.
The different approaches of the two collaborations means that they are not measuring the same quantity so that the values obtained cannot be directly compared to one another. The predictions for the measurements made by each collaboration for a given hadronic model and composition are also different. In spite of this, the two sets of measurements [43,44] have been superposed in the right panel of Fig. 13. This figure must be interpreted with great care. The two sets of measurements are compared to predictions for protons for the true < X max > using EPOS-LHC [45], SYBILL2.1 [46] and QGSJETII-04 [47] (the three upper lines) and for the effective < X max > expected in the Telescope Array, using QGSJET-01 [48], SYBILL2.1 and QGSJETII-03 both for proton and iron primaries (respectively the middle and lower sets of lines). We note that the most recent hadronic models EPOS-LHC and QGSJETII-04, and SYBILL2.1 have systematically deeper X max than earlier versions of QGSJET and EPOS. Predictions for protons for the two experiments with SYBILL2.1 illustrate the magnitude of the selection bias of the Telescope Array data. Clearly the Pierre Auger measurements of X max are compatible with 18pc Figure 13. Effective < X max > measured at the Telescope Array (small black squares) and < X max > obtained from the Auger Observatory data with fiducial cuts (large brown circles) compared to expectations for the effective X max measured by the TA for proton (middle set of lines; blue in on-line version) and iron (lower set of lines; red in on-line version), using QGSJET-01 (dot dashed), QGSJETII-03 (full) and SYBILL-2.1 (dashed). The upper set of lines (purple in on-line version) correspond to X max for proton with EPOS-LHC (dot dashed), SYBILL-2.1 (dashed) and QGSJETII-04 (dot long dashes) from top to bottom (see text).
a proton like composition at energies of order 1 EeV. The rate at which the shower maximum changes with energy, or elongation rate, is much smaller than expectations for a pure composition suggesting that composition is becoming heavier. The Telescope Array results are, in contrast, compatible with expectations for protons using QGSJETII-03.

Combining the average and variance of the depth of maximum
Study of the X max distributions has not helped to clarify this controversy. The variance measured by the Pierre Auger Collaboration has consistently shown a decrease with energy above 1.8 EeV which is suggestive of a composition becoming both heavier and less mixed. The TA results indicate distributions in effective X max which are compatible with proton simulations. It is worth noting that the total number of points in the TA analysis with energy exceeding 1.6 EeV (438) is of the same order as those reported in the single bin between 5 and 6.3 EeV (413) in the last update by the Pierre Auger Collaboration [49]. The implications of these measurements depend on the interaction model. The mean value of the X max distribution is indicative of the average mass composition and the variance can be related to the dispersion in composition. This can be exemplified using a simple generalization of the Heitler model for shower development [50]. In this simple model it can be shown that the average and variance of X max are given by: Here the sub-index p refers to a pure proton composition and f E is an energy dependent term. The variance of X max has a contribution from the variance of each primary, σ 2 sh , and another because each primary has a different < X max >. Here σ sh can be regarded as a function of mass number, A, and can be parameterized as: Taking the average of Eq. 3 it is straightforward to invert the above equations to give an explicit relation of the average and variance of lnA in terms of measurable quantities: These relations depend on the hadronic model though < X max > p , f E and < σ sh >.
The results for three different models are shown in Fig. 14. The averages obtained of lnA indicate a clear transition from light to heavier composition as the energy rises above 2 EeV. The variance of lnA indicates the level of mixing and a value of zero corresponds to a pure composition no matter what. The negative variances obtained for QGSJETII-04 are unphysical, suggesting that either the model is not compatible with the data or that the measurements are unreliable [50]. These results have been recently updated and completed with an analysis of the full X max distributions, which allows estimates of the relative abundances assuming a mixed composition of proton, helium, nitrogen and iron [49]. The conclusions are basically the same.

The distributions of depth of maximum
The analysis of the mean and variance does not fully exploit the quality and richness of the available data which can be obtained from the distributions of X max . The Telescope Array has compared the effective X max distributions obtained in different energy bins and has found them compatible with proton composition and QGSJETII-03 [40]. The distributions obtained by the Pierre Auger Collaboration have been fitted to obtain the fractions of primaries according to different hadronic models. Templates are obtained for distributions of a restricted number of pure primaries, namely proton, helium, nitrogen and iron. Best fit values of the relative fractions of the primary composition have been obtained in three simplified scenarios, namely restricting the number of different species to proton and iron only, considering proton, helium and iron and finally using a four component mix by also including nitrogen [41].
The first scenario with only proton and iron is shown to give poor quality fits but the situation improves as intermediate mass primaries are introduced. The results obtained with EPOS-LHC are satisfactory for all the energy range while QGSJETII-04 does not fit the distributions well. It can be argued that adding extra components will not improve the fit for this model. This is consistent with the unphysical results obtained for QGSJETII-04 when the first two moments (mean and variance) of the X max distribution are interpreted in terms of the mean and variance of lnA as discussed earlier. The predictions of proton and iron fractions are remarkably similar   Figure 14. Results for < ln A > (top row) and σ 2 ln A (bottom row) inferred with the Auger X max data for three different hadronic models.
for all models while they differ greatly on the composition of intermediate mass ranges. The proton fraction is energy dependent changing from values between 60% and 80% at energies around 2 EeV but drops to very low fractions as the energy increases to 10 EeV. The proton fraction in the EeV range seems to be high to explain the ankle feature with the transition to extragalactic sources unless there is a proton contribution at EeV energies from the galaxy too. But this is disfavored by anisotropy measurements as explained in the previous section. On the other hand a pure proton scenario would require substantial changes to the hadronic models.

Measurements of the muon content of showers
Highly inclined showers consist mainly of muons and thus provide an opportunity to measure the muon content of air showers [51]. Since N 19 is the relative size with respect to a 10 EeV proton shower simulated with QGSJETII-03, it is also an approximate estimate of the relative number of muons with respect to the same reference. Detailed comparisons have been made between the number of muons relative to those of the reference and the fitted values of N 19 in simulated showers using different composition and hadronic models. They reveal a small energy dependent bias and spread. A new variable R µ is defined correcting N 19 at the level of a few percent to compensate this bias and to minimize the dispersion from hadronic models [52]. The calibration data for inclined showers provides a measurement of R µ and an independent measurement of the FD energy. The results are also consistent with a change to heavier composition in the EeV range [51]. There are more muons in the data that in the reference by about a factor of 1.8. This can be marginally accommodated using iron and the most recent hadronic models. However measurements of X max favor a much lighter composition which is in tension with these results. This is illustrated in Fig. 15 that displays a two-dimensional plot of < X max > and ln R µ at an energy of 10 EeV compared to expectations from different hadronic models and mass compositions.
Similar results have been observed when comparing the energy of events with θ < 60 • obtained  Figure 15. Two dimensional plot of ln R µ obtained from inclined showers and X max at the Pierre Auger Observatory compared to expectations form different models and different primaries. Figure 16. Measurements of the average value of X µ max as a function of energy using the Pierre Auger Observatory SD data. The results are compared to expectations for both proton and iron using the most recent hadronic models.
using simulations to those obtained with the calibration procedure and with other methods to count muons [53,54]. This is not new, similar evidence has been found from experiments measuring horizontal muon bundles [55], from lower energy surface arrays such as KASCADE-Grande [57] and from the Large Hadron Collider (LHC) [56].

Measurements of the muon production depth
The hybrid approach has great potential to explore new SD observables in order to further constrain composition and hadronic models. The Pierre Auger Collaboration has explored observables related to the muons in the shower front because it is very sensitive to them. An interesting method has been developed by studying the time structure of the muon front [58]. The time delays of the muons with respect to a plane front can be shown to be dominated by geometry and a simple relation holds between the lateral displacement, the time delay and the distance at which the muon was produced. Under some conditions, this relation can be inverted to obtain the distribution of the depths at which the shower muons are produced, using the time structure of the data of the SD stations [6]. This variable has been designated as X µ . The depth at which the maximum production is achieved, X µ max , depends on composition and hadronic models as could be expected, similarly to X max .
The average values of X µ max measured by the Pierre Auger Observatory are displayed as a function of energy in Fig. 16 compared to expectations from several models. Interestingly, one of the required conditions to obtain it is that the number of sampled muons is high, so that the method is more sensitive at the highest energies where the hybrid data are scarce. The results are also consistent with the primary mass composition becoming heavier. However the values obtained are out of the bracketed values obtained with EPOS-LHC. The measured value of < X µ max > can be also converted to an average value of < lnA > mimicking the procedure used for < X max >. The results for the values of < lnA > obtained using both X max and X µ max are combined in Fig. 17 for QGSJETII-04 and EPOS-LHC. The two most recent hadronic models QGSJETII-04 and EPOS-LHC, updated using LHC data, seem to be inconsistent with X max  Figure 17. Results for < ln A > obtained using the measurements of both X max and X µ max data at the Pierre Auger Observatory using the most recent hadronic models. and X µ max data taken together, pointing to the possibility that part of the puzzle lies in the extrapolations of accelerator data. Surely this is one of the paths that must be further explored in the future [59].

Summary and outlook
Recent results obtained by the TA and the Pierre Auger Collaborations have been reviewed. The spectrum measurements agree below 50 EeV and display the ankle at about 5 EeV. For energies greater than 50 EeV, both spectra clearly indicate a suppression but there are some important discrepancies. TA results are consistent with models in which protons are accelerated and they interact with the CMB producing the GZK cutoff. The Auger spectrum can only be made consistent with such scenario if the fluxes at the sources are suppressed at similar energies. But such a scenario is in fact in tension with X max and X µ max data from the Auger Observatory. Anisotropy searches with the Auger data have shown some evidence of large-scale anisotropy which is becoming stronger and both observatories have indications of possible hot spots in each hemisphere. The significance of the TA hot spot is stronger. Recent updates of these measurements enhancing the sky coverage with inclined showers and combining the data of the two observatories have given compatible results strengthening the case for large-scale anisotropies. It is in fact not ruled out that large scale anisotropies may have some impact in the flux discrepancies between the northern and southern hemisphere.
Composition measurements also indicate discrepancies between the two experiments. Unbiased X max measurements can be related to primary mass if particle interactions are under control. Unfortunately this is not the case and hadronic models are evolving as new data appear. While the TA group measures an effective X max because of selection effects, the Pierre Auger Collaboration has developed fiducial cuts to avoid bias. Results are not easy to compare since each group uses observables not strictly identical and compares to predictions using different hadronic models. The TA flux and X max data are claimed to be consistent with a pure proton scenario. The average and variance of X max Auger data indicate a transition from light to heavier elements starting at E ∼ 4 EeV. We note that at roughly this same energy the spectrum displays the ankle and there are indications of a small anisotropy signal. Auger data allow the measurement of X µ max and the numbers of muons in inclined showers (and with other methods). These data also suggest that the primary mass composition becomes heavier as the energy increases, however new difficulties arise when all the data are tried to be interpreted in terms of conventional hadronic models. The measured number of muons in the showers is difficult to reconcile with the deduced composition inferred from X max measurements. On the other hand the interpretation of X µ max and X max results in terms of composition is inconsistent for EPOS-LHC. The average and variance of the X max distributions measured at the Auger Observatory disfavor QGSJETII-04, giving a negative variance of lnA. Older models do not improve the picture. Finally other interesting findings can be attributed to the two groups, such as negative searches for photons, neutrons and neutrinos that have led to relevant bounds.
Auger results related to mass composition using these independent methods suggest a mixed composition in which, qualitatively, nuclei of increasing A are dominating the composition as the energy increases. When these results are considered together with the energy spectrum, new classes of models have been proposed such as those supporting a heavy composition at the sources which gets degraded through interactions with the CMB, or those which assume that the sources have a maximum acceleration energy which increases with the nuclear charge, Z [20,60,61,62,63]. But these interpretations must be taken with caution because standard hadronic models fail to interpret all the data in a consistent way.
While much progress has been made in the study of UHECR there are still many fundamental issues to be solved. The different interpretations obtained by the Telescope Array and the Pierre Auger collaborations should be understood. There are several joint groups working together on the comparison of the data and combined analyses, notably for spectrum, composition and anisotropy studies. There are plans for new enhancements and extensions of the two observatories, to increase the size of the TA [64] and to add muon counters at the Auger Observatory. The results obtained with these enhanced detectors and the collaborative efforts of the two collaborations should clarify the picture in the next few years of what is likely to shape the future of the field.