Resonant Diphoton Phenomenology Simplified

A framework is proposed to describe resonant diphoton phenomenology at hadron colliders in full generality. It can be employed for a comprehensive model-independent interpretation of the experimental data. Within the general framework, few benchmark scenarios are defined as representative of the various phenomenological options and/or of motivated new physics scenarios. Their usage is illustrated by performing a characterization of the 750 GeV excess, based on a recast of available experimental results. We also perform an assessment of which properties of the resonance could be inferred, after discovery, by a careful experimental study of the diphoton distributions. These include the spin J of the new particle and its dominant production mode. Partial information on its CP-parity can also be obtained, but only for $J\geq2$. The complete determination of the resonance CP properties requires studying the pattern of the initial state radiation that accompanies the resonant diphoton production.


Introduction
The resonant production of a photon pair at hadron colliders is quite a simple process, which we can hope to characterize with a high degree of generality. To do so, first of all we need to understand the possible initial states that can lead to the production of the intermediate resonance R decaying to γγ. If no additional hard objects are present in the final state, which is our working hypothesis, only a few partonic scattering processes are likely to be relevant, namely the ones involving gluons (gg), quarks (qq, with q = u, d, c, s, b) or photons (γγ). "Mixed" situations such as qg-initiated production are forbidden by color conservation and by Lorentz symmetry, which requires the heavy resonance R to have integer spin J (with J = 1 by the Landau-Yang theorem). Channels of the type q q with q = q are strongly disfavored by flavor constraints, which make very difficult to imagine how a resonance within the energy reach of the LHC might have sizable flavor non-diagonal couplings to the light quarks. We will thus ignore this possibility in what follows.
Among the other channels that might be considered we can definitely exclude tt, because tt-initiated production unavoidably comes together with a tt pair in the final state from the splitting of the initial gluons, while we choose to limit our analysis to final states with no extra hard objects. Although similar considerations hold for bb production, the associated b-quarks are typically soft. Thus they are hard to detect and to identify as b-jets and can be easily confused with the radiation pattern that characterizes the other partonic modes. Still, after the identification of a signal and with large enough luminosity, checking for the presence of bottom quarks will allow to distinguish the bb mode from the others.
Production through Massive Vector Bosons (MVB), namely W + W − -, ZZ-or γZ-initiated processes, will also be neglected. 1 The MVB processes are marginal to the present study for two reasons. First, they are accompanied by the production of forward energetic jets from the quark splitting, which are typically hard enough and not too forward to be detected. MVB production is thus distinguishable from the partonic processes provided suitable forward jet selection cuts are put in place. Notice that the situation is different for the γγ production mode because the photon is massless and thus the p ⊥ of the emission is only cut-off by the proton mass. The QCD jets from γγ fusion are thus softer than the MVB ones and difficult to detect. Actually, the γγ radiation pattern is even softer (and possibly even consist, in the elastic scattering regime, of just two extremely forward protons) than the one associated with the other partonic processes gg and qq, giving a possible handle to pin it down [2]. The second reason to neglect the MVB processes is the fact that the photon parton distribution function (PDF), again because of the lack of a hard low-p ⊥ cutoff in the photon splitting, is larger than the MVB one. 2 This makes MVB processes also quantitatively marginal. An exception is the situation in which the couplings of R to MVB's are much larger than the γγ one, in which case, however, resonance searches in MVB final states are much more effective.
On top of the analysis of the possible production channels, a full study of a resonant diphoton process also requires a characterization of the cross section and kinematical dis- The incoming partons (of helicity λ 1,2 ) annihilate into a resonance of spin J (and spin-projection m = λ 1 − λ 2 along the beam axis) that subsequently decays into two photons with helicities λ and λ . We denote as S = |λ − λ | the absolute value of the spin of the diphoton system along the decay axis.
calculations of the signal rate and distributions in terms of the parameters that control the on-shell resonance production and decay Feynman amplitudes. The translation of the latter parameters into effective operator coefficients, which straightforwardly allows to implement our signal in an event generator in order to deal with QCD radiation and detector effects, is reported in Appendix A for J = 0 and J = 2 resonances. In the fully general case, in which all the 7 production modes are active and no further assumption is made on the resonance couplings, the proliferation of free parameters makes the problem untreatable. Therefore in section 3 we define a set of representative benchmark scenarios, whose number and variety should be sufficient to provide a wide enough coverage of the various phenomenological options. These scenarios are analyzed by recasting, with a strategy described in Appendix B, available 8 and 13 TeV experimental searches. Rather than aiming to fully quantitative results, which might be only obtained by the experimental collaborations, the goal of this study is to illustrate the usage of our benchmarks to characterize possible signals such as the popular 750 GeV excesses. Still, we will be able to reach semiquantitative conclusions on the viability of our scenarios. In section 4 we report our conclusions and a preliminary assessment of the additional information which can be extracted from the study of initial state radiation emission. A complete discussion of the latter point is left for future work.

General framework
We consider a resonance R of integer spin J, produced at the LHC out of a given 2-partons initial state in = {gg, qq, qq, γγ} and decaying to γγ as depicted in figure 1. 4 We start our discussion from the fully polarized scattering process and denote by λ 1 and λ 2 the helicities of the incoming partons, λ and λ those of the final state photons. These helicities cannot be measured at the LHC and we will eventually have to sum/average over them to obtain the cross-section. 5 Conservation of angular momentum along the beam direction implies 4 Initial partons are ordered by the direction they come from, this is why qq and qq are distinct in states. 5 We assume that it will never be possible to measure photon polarizations at the LHC and we restrict our attention to inclusive γγ production. The exclusive case, in which we imagine having access to the radiation that only a single spin component of the resonance can contribute to the partonic process, namely the one with spin projection m = λ 1 − λ 2 along the beam axis oriented in the direction of parton "1". Thus, the resonance production process can be fully described by a set of dimensionless coefficients A in λ 1 λ 2 which parametrize the corresponding Feynman amplitudes as where M is the resonance mass. The helicities λ 1,2 can assume the values λ 1,2 = ±1 for in = gg or in = γγ and λ 1,2 = ±1/2 for in = qq/qq. Correspondingly, only resonances with m = 0, ±2 and m = 0, ±1 can be produced, respectively, in the bosonic and fermionic channels. Not all the four complex coefficients A in λ 1 λ 2 we have for each production mode are independent. Invariance of the amplitudes under a π rotation in a direction orthogonal to the beam implies A gg/γγ where in the first equality we implicitly made use of the fact that the gg and γγ states are made of indistinguishable particles.
In the case of the resonance production, which we discussed until now, the incoming partons momenta are completely fixed in the COM frame, thus it is trivial that the amplitudes can be parametrized in terms of few constant coefficients. The situation is different for the resonance decay process, which depends on the kinematical variables of the γγ final state and in particular on the COM scattering angle θ. Still, each polarized decay amplitude can be parametrized by a single constant because the angular dependence is completely determined, and encapsulated in the so-called "Wigner d-matrices" d J m,m (θ) [11]. The point is that by a rotation one can connect a γγ state with arbitrary polar and azimuthal angles θ and φ to a photon pair moving along the beam axis and obtain the angular dependence of the amplitude from the matrix elements of the rotation matrix among the resonance spin eigenstates. The result reads (see for instance [7]) where we made use of the CPT symmetry to relate (up to phases, which eventually produce the (−) J factor) the amplitude coefficients of the R → γγ decay to those associated with the production process γγ → R. Therefore describing the resonance decay does not require introducing new parameters. The set of processes we are considering is thus fully characterized, taking into account the relations in eq. (2), by a rather small number of parameters shown in table 1. Namely, we have in general 4 complex parameters for the qq (and qq) production, 3 complex parameters describing gg/γγ if J is even and only 1 complex parameter if J is odd. For J = 0, the "+−" and "−+" amplitudes vanish and we are left with 2 complex parameters for gg/γγ and again 2 for the qq channels. The case J = 1 is not worth discussing because the decay to γγ (and the production from gg) is forbidden by the Landau-Yang theorem, or equivalently by noting that also a g/γ 2 vanishes in this case (see table 1) because of angular momentum conservation.
It is important to remark that the derivations above are completely model-independent as they only rely on the invariance under rotations and CPT, which are symmetries of any from the initial state, is briefly discussed in section 4. Table 1: Amplitude coefficients expressed in terms of a set of complex parameters "a". Untilded and tilded parameters are, respectively, CP-even and CP-odd. For shortness, +1 and +1/2 helicities (which are appropriate for gg/γγ and qq initial states, respectively) are both denoted as "+" and the same for "−".
relativistic quantum theory of particles. In particular they do not rely on the CP symmetry, therefore our results hold irregardless of the resonance CP-parity and even of whether CP is at all a symmetry or not. If CP is a symmetry, we get the additional constraint where ρ CP = ±1 is the intrinsic CP-parity of the resonance. Therefore only some of the parameters, denoted as untilded a's in table 1, survive for a CP-even resonance and only tilded ones in the CP-odd case. Sizable tilded and untilded parameters would be simultaneously present only if the CP symmetry was badly broken by the resonance couplings. We stress that the "a" (and " a") coefficients in table 1 are, in general, complex numbers. 6 However they become real when the resonance production/decay processes are induced by heavy mediators. Establishing experimentally whether they are real or not would therefore allow us to verify or falsify this hypothesis. In order to appreciate this claim, we notice that if the resonance couplings are mediated by the exchange of heavy particles it is possible to integrate them out, giving rise to a set of local operators (contact interactions) that induce resonance production and decay. The heavy-mediator condition can thus be equivalently formulated as the hypothesis that the production/decay amplitudes are well described by a contact interaction at Born level, i.e. by the matrix element of a local Hermitian operator, in which case the CPT symmetry, combined with eq. (2), gives a relation It is easy to check that this condition implies that the a's in table 1 are real. If instead the resonance couplings are due to light particles loops, imaginary parts will arise in the amplitudes, by the optical theorem, due to the propagation of on-shell intermediate states.
Establishing whether the a's are real or not would thus give us relevant information on the 6 In spite of the fact that they were erroneously taken real in the first version of the manuscript. We thank R. Rattazzi for pointing this out to us. resonance dynamics. However, this will turn out to be impossible through the measurement of unpolarized inclusive diphoton production distributions, which are our main target. Indeed, restricting to real a's is enough to span the whole variety of kinematical distributions one would obtain even for general complex a's. We will briefly come back to this point in section 4.
A priori, the parametrization of the resonance production and decay amplitudes provided in table 1 might still be redundant because they solely followed from rotation and CPT invariance. In principle, further constraints might arise by requiring invariance of the amplitudes under the complete Lorentz group. This is however not the case, as explicitly shown in Appendix A for J = 0 and J = 2 resonances. 7 In the Appendix, we classify all the Lorentz-invariant terms, expressed as functions of the 4-momenta of the resonance and in particles and of their polarization vectors or spinor wave functions, which can appear in the polarized amplitudes. The coefficients of these Lorentz-invariant terms are found to be in one-to-one correspondence with the parameters in table 1, showing that no further restrictions emerge from imposing the full Lorentz symmetry. Moreover, Lorentz-invariant amplitudes are easily mapped to Lorentz-and gauge-invariant operators and therefore another result of the Appendix (see eq.s (31) and (32)) is to relate the phenomenological parameters a i , a i , and in turn the A in λ 1 ,λ 2 's, to the couplings of a phenomenological Lagrangian. 8 This is required for the implementation of our parametrization in a multi-purpose event generator. Consistently with the discussion following (5), if the phenomenological Lagrangian is taken to be Hermitian (i.e. the only phenomenologically relevant states are R, in, γγ) the amplitude coefficients obey eq. (5) and the a's are real, as expected. In the Appendix we focused on J = 0 and J = 2 resonances because higher spin particles are anyhow not implemented in multi-purpose event generators. Complete simulations for J ≥ 3, taking properly into account soft QCD radiation, hadronization and detector effect would thus require a different approach, based on matrix-element reweighting techniques as discussed in the next section.

Partonic cross sections
We are now in the position of constructing, with the amplitude coefficients as building blocks, the partonic unpolarized cross-section of the complete 2 → 2 reaction in → γγ. This will allow us to identify the combinations of amplitude coefficients that appear in the unpolarized cross-section and will suggest a convenient phenomenological parametrization of the signal, to be employed for the experimental characterization of the resonance properties.
The in → γγ Feynman amplitude is the product of the production and decay amplitudes, times the Breit-Wigner propagator of the resonance where Γ is the total resonance width andŝ is the partonic COM energy squared. The Breit-Wigner propagator produces, in the amplitude squared, a factor of π/(ΓM ) times the 7 The case J = 3 has also been checked, but it is not discussed in the Appendix. 8 Notice that the correspondence among the Lorentz-invariant terms in the (on-shell) amplitude decomposition and the operators is not at all one-to-one. Namely, infinitely many operators reduce, on-shell, to a single term in the amplitude. The simplest set of operators, just sufficient to produce arbitrary on-shell amplitudes, is selected in the Appendix.
normalized Breit-Wigner distribution BW(ŝ). The partonic cross-section thus reads having reabsorbed inσ in some factors, and in particular the dependence on Γ. The second equality in the equation holds for a narrow resonance, namely in the limit Γ/M → 0. In that limitσ in assumes, as we will readily see, the physical meaning of the signal cross-section for unit parton luminosity at the resonance mass, namely for [τ dL in /dτ ]| τ =M 2 /s = 1. A compact expression for dσ in /d cos θ (see eq. (9) below) may be obtained as follows.
As previously explained, each polarized in → γγ process is mediated by a single resonance spin m = λ 1 − λ 2 . Therefore its angular dependence is fixed by the Wigner formula (3) to be the square of the associated Wigner matrix, [d J m,λ−λ ] 2 . By summing the polarized cross sections over m, λ and λ we obtain the unpolarized one, expressed as a sum of known functions of the COM scattering angle θ weighted by the square of the corresponding polarized production and decay amplitudes. The polarized production amplitudes can be traded for the resonance production cross section, whereas the decay amplitudes can be traded for the branching ratios.
In the sum, several terms can be grouped together by proceeding as follows. We first sum over the photons helicities λ and λ and notice that the ++ and −− terms in the sum produce the same angular function, The functions D (J) |m|,S (θ) have also appeared in previous work, see e.g. [10]. Here we chose to normalize them to unity in the integration domain cos θ ∈ [0, 1], which is the appropriate one since the final state photons are indistinguishable particles. The above equation ensures that terms in the λ 1,2 sum with a given value of m = λ 1 − λ 2 have the same angular dependence of those with the opposite value, so that the two can be grouped in a single term. The double sum over the initial state polarization thus becomes a single sum over the absolute value of m, |m| = 0, 1, 2.
Since S = 0, 2 ranges over two values and |m| = 0, 1, 2, six terms are present in the sum, each characterized by its own angular distribution D (J) |m|,S (θ). Notice however that only four of the six terms can be simultaneously turned on in a given partonic process because |m| = 0, 1 for in = qq and |m| = 0, 2 for in = gg/γγ. Nevertheless we will momentarily retain the six of them for a more concise exposition.
The unpolarized cross section can finally be written as  0,2 , leading to only five independent distributions. Moreover, since the only viable values of |m| are 0, 2 for gg/γγ production and 0, 1 for qq, only three distributions are present in the former case and four in the latter. 9 For J = 2, the distributions relevant for gg/γγ and for qq are displayed in the plots in figure 2. We see they have considerably different shapes so that it should be possible to distinguish them even with moderate experimental accuracy.
The cross sections and branching ratios appearing in eq. (9) are defined as follows. Thē σ's are the total cross sections (at unit parton luminosity) for the production of the resonance with spin m = |m| plus, if m = 0, the one with m = −|m|. Namelȳ where c in = 1, 3, 8 are the color factors for, respectively, in = γγ, qq, gg and |A gg/γγ 2 The cross-section for in = qq need not to be discussed explicitly because it is just identical to the qq one by the second relation in eq. (2). The BR's in eq. (9) are those for the resonance decaying to a polarized diphoton pair with equal helicities λ = λ = ±1 for S = 0 and with opposite helicities for S = 2, i.e. with A γγ S again as defined in eq. (11). Notice that the fact of having two distinct decay channels (++ and −−) for S = 0 and only one (+−, which is indistinguishable from −+ after angular integration) for S = 2 compensates for the fact that the ±± states are made of indistinguishable particles and thus they have to be integrated over half of the solid angle. Furthermore, the branching ratios, as apparent from the notation, do not depend on the resonance spin m because of rotational invariance.
Eq.s (10), (11) and (12) provide the required map among the amplitude coefficients and the potentially observable quantities (σ's and BR's) that parametrize the partonic cross section in eq. (9). We see that the observables depend on few combinations of the a and a parameters that control the amplitude coefficients through table 1. In particular, no information can be extracted on whether the a's are real or complex, namely on whether eq. (5) is satisfied or not, as previously mentioned.
The cross section parametrization in eq. (9) can be directly employed for the comparison with experiments or, as we find convenient to do for the analysis in section 3, rewritten in a "probabilistic" format by factoring out the total resonance production cross section times the total branching ratio to an unpolarized photon pair (BR), namely Here is the probability for the produced resonance to have spin equal to |m| in absolute value and to decay to a state of spin S. The last identity in eq. (14) has been obtained using eq.s (10) and (12). The probabilistic format is useful as it allows to disentangle the total signal rate from the normalized angular distribution, encapsulated in the P's. Notice that the P's, precisely because they are probabilities, sum up to one.

LHC cross sections and distributions
It is conceptually straightforward to go from the partonic cross section, characterized by the σ×BR and P parameters as in eq. (13) (or by eq. (9)), to LHC differential cross sections or to event samples to be compared with the experimental data. The result will consist in a linear combination of distributions or in an admixture of event samples, each generated with its own "D" partonic distribution and weighted by the corresponding "P" probability. Such event samples could be obtained in two ways. Either by direct simulations, from the Lagrangian in Appendix A implemented in MadGraph [12], turning on at each time the couplings associated with a given "D", or by matrix element reweighting, starting from the simulation of a scalar and reweighting each event, with partonic scattering angle θ, by D(θ). This latter approach is the only viable one for J > 2, where no multi-purpose event generator implementation is available as previously mentioned.
For an accurate comparison with the data, properly taking into account soft QCD radiation, hadronization and detector effects one of the two strategies described above should be adopted. For the illustrative purpose of the present paper, however, it is sufficient to stick to purely leading order predictions, on top of which experimental effects will be attached as overall efficiency factors as described in the next section. This simple approach has the advantage of producing semi-analytical formulas for the distributions from which we can get an idea of which aspects of the signal properties are easier to extract from data.
The cross section, differential in the cosine of the scattering angle in the COM and in the boost of the COM frame, reads having made use of the right hand side of eq. (7), that holds in the narrow resonance limit Γ/M → 0. In the above equation, τ = M 2 /s, with "s" the collider energy squared, τ dL in /dτ is the differential parton luminosity and dP in /dy is the distribution of the COM boost y. These functions are related to the initial state PDF's f by The variables y and cos θ are related to the rapidity of the two photons and to their p ⊥ as  Figure 3: Differential parton luminosities dP in /dy, as defined in eq. (16). The left plot shows gg, uu and dd initial states while cc, bb, ss, γγ (and again gg for comparison) are displayed on the right. The 1 σ bands are obtained as described in the text.
Notice that cos θ ranges from 0 to 1 as the photons are indistinguishable.
Both cos θ and y are measurable and both the cos θ and y differential distributions contain interesting and, to large extent, complementary information about the signal. Namely, the cos θ distribution gives us direct access, at least if only one "in" channel is active in eq. (15), to the partonic differential cross section, which in turn is related to the resonance spin as previously discussed. It also provides partial information about the production mode, given that the cos θ distributions, i.e. the D functions, can be different if the resonance is produced by the gg/γγ or by the qq initial state. 10 It is instead unable to distinguish, for instance, qq = uu from qq = dd as the D's are the same in the two cases. The situation is basically reversed for the differential distribution in y, which is insensitive to the details of the partonic cross section and is entirely dictated by the production mode, which determines the shape of dP/dy. Whether or not and how easily this may be exploited to distinguish different production mechanisms depends on how much different the dP/dy's are in the different cases. This is quantified in figure 3, for a resonance mass of M = 750 GeV (chosen in preparation for the discussion of the next section) and √ s = 13 TeV. We see that the two valence quarks have slightly different distributions, allowing in principle to distinguish uu from dd. All the sea quark distributions are instead very similar, or even identical within the uncertainties, and not far from the ones for gg and γγ. The plots in figure 3 are obtained by the NNPDF23 nnlo as 0119 qed set of NNPDF2.3 [4] with a factorization scale of 750 GeV. The uncertainties are obtained from the variance over the PDF replicas provided in the PDF set. Scale uncertainties, quantified by varying the factorization scale, are found to be negligible. This is valid for the "ordinary" partons g and q, but not for the photon, whose PDF measurement is too bad to extract any quantitative information. The γγ luminosity is thus taken from ref. [2], where it has been estimated from the theoretical calculation of the photon PDF presented in refs. [3,13]. Uncertainties in dP γγ /dy are not reported in ref. [2] and consequently they do not appear in our plot. 10 However they can also be equal, since we saw in the previous section that D (2) 0,0 and D (2) 0,2 can appear in both gg/γγ and in qq production. If this is the case, distinguishing the two channels requires looking at the y distribution as we will readily discuss.  We saw that cos θ and y differential distributions provide complementary information about the signal, however because of the photon acceptance cuts it is not clear that the two distributions can actually be disentangled experimentally and measured separately. While performing separate measurements (possibly unfolding the experimental effects) would facilitate the interpretation, allowing for instance to compare directly the cos θ distribution with the shape of the D functions in figure 2, notice that the exact same amount of information could be extracted from the study of the doubly differential distribution.
Before concluding this section and in preparation for the next one, where we will use our framework for a first characterization of the 750 GeV excess, we report in table 3 the total parton luminosity at M = 750 GeV at the 13 and 8 TeV LHC and the gain, defined as the ratio of the 13 and 8 TeV cross sections, for each production mode. Contrary to dP/dy, the uncertainties are now dominated by scale variation and is of order 10% (up to 15% for the gluon, and below 6% for quarks). Two set of results are reported in the table concerning the γγ channel. The first one is based on the theoretical prediction from Ref. [2], which we will employ in what follows. The second one, subject to a large error, is obtained with the NNPDF2.3 [4] PDF set.

Benchmark scenarios
In the previous section we saw how the production of a resonance of arbitrary spin decaying to γγ is conveniently parametrized, for each given in = gg/γγ/qq production channel, in terms of a rather small number of phenomenological parameters with a sharply defined and intuitive physical meaning. However, being completely agnostic about the resonance couplings would require taking all the production channels into account simultaneously, with independent free parameters for each of the 7 (i.e., gg/γγ/qq = {uu, dd, cc, ss, bb}) in states. This proliferation of parameters makes the problem untreatable in full generality and obliges us to make additional assumptions in order to reduce the dimensionality of the parameter space. A set of plausible restrictive assumptions is defined in the present section, producing a set of alternative benchmark scenarios. Each of these benchmarks contains a small enough number of free parameters to be experimentally tested in full generality. The variety of benchmarks should provide a sufficient (but still unavoidably partial) coverage of the phenomenology. Additional benchmarks can be defined, if needed, within our general framework.
The benchmark scenarios can be used for exclusions, producing limits on σ×BR which are more general and easier to reinterpret in specific models than those obtainable with the habitual benchmarks of a scalar or of a J = 2 "RS graviton" resonance. More interestingly, they can be used to characterize the properties of a new resonance that we might happen to discover in the diphoton final state. In the latter case, the SM p-value and other statistical quantities aimed at assessing the actual existence and viability of the signal, could be reported on the benchmark model parameter space. This will select the signal hypothesis that best fits the data and will give us information about the resonance spin and (see section 4) CP properties. At a later stage, with enough data, it will be possible to measure the parameters of the benchmark models, namely those that control the signal kinematical distribution and the total σ×BR. The model-independent nature of our parametrization will straightforwardly allow to translate these measurements into whatever the "true" resonance model turns out to be.
A good fraction of the program outlined above is slightly premature, as a discovery still has to come. However the M = 750 GeV excess reported by ATLAS [14] and CMS [15] with 13 TeV LHC data gives us the opportunity to practice, at least on some aspects of the signal characterization strategy. 11 We will do so by recasting ATLAS 13 TeV [14], CMS 13 TeV [15], ATLAS 8 TeV [18] and CMS 8 TeV [19] experimental searches, with a procedure described in Appendix B in detail. It suffices here to say that the recast is performed by reconstructing, in the Gaussian approximation, the likelihoods associated to each experimental search from the background-only p-value and the observed limit. The four searches are eventually treated as statistically independent in the combination. The intrinsic inaccuracy of our statistical method and our approximate treatment of the experimental efficiencies make our results not fully quantitative. Moreover, the experimental searches we use are not optimized to provide information about the angular distributions of the putative signal and thus they are poorly sensitive to the resonance spin and production mode. Consequently our results will often show a rather limited discriminating power within the parameter space of each benchmark and among different benchmarks. Most of what we will be able to tell will come from the combination of 8 and 13 TeV searches because of the slightly different gain factors r (see table 3) in the total signal rate. Notice however that the situation would substantially improve with dedicated experimental analysis and/or more data.
In view of the considerations above, we warn the reader that the results that follow should be mostly regarded as a pragmatic illustration of the usage of our benchmarks. Still, it will be interesting to see that in some cases the various analyses do display slightly different acceptances for the same signal shape, merely due to the slightly different selection cuts. This produces, in the combination, some discriminating power among the different hypotheses and indicates that progress in the signal characterization should be relatively easy to achieve with a dedicated analysis.

Scalar resonance
As a first case we consider the simplest scenario, that is the model with a scalar resonance. This case is rather peculiar since the angular distribution of the two photons in the COM frame is completely flat. Indeed, as we saw in the previous section, the only contribution to the production comes from the m = S = 0 mode, which is described by the angular function D (0) 0,0 = 1 (see table 2). The only model-dependence is encoded in the relative strengths of the 11 Provided that the signal originates from a single resonance decaying in a photon pair rather than a pair of axions decaying into highly collimated photons as suggested in ref. [16] (see also refs. [17]). various production channels, which can be parametrized through the partonic production cross sectionsσ in . Such a parametrization characterizes the possible scenarios in a way that is completely independent of the details of the experimental searches, in particular of the COM energy of the collider. From a practical point of view, however, this does not seem a convenient choice. Due to the extremely different parton luminosities (see table 3), partonic cross sections of similar size give rise to signal cross sections for the various production channels that can differ by more than one order of magnitude. For instance, the production modes through quarks or photons can be comparable to the gg one only if their partonic cross sectionsσ qq/γγ are much larger thanσ gg . Therefore, we find more convenient to adopt a parametrization that allows to efficiently treat cases in which various production modes lead to comparable signal yields. Of course a parametrization of this kind is necessarily collider dependent, since it must take into account the parton luminosities. A possible choice, which we will adopt in the following, is to use the ratios of signal cross sections for the various channels for a collider energy of 13 TeV. In particular we define the quantities where σ 13 TeV in is the 13 TeV production cross section in the in channel, whereas σ 13 TeV tot is the total production cross section. 12 The R in parameters directly encode the relative importance of the contributions to the signal cross section from the various production channels. Since they are normalized to the total production cross section, the R in parameters sum up to unity, in R in = 1. The relative strengths of the production channels at 8 TeV can be easily related to the 13 TeV ones by taking into account the change in the partonic cross sections listed in table 3.
From the experimental point of view, the various production channels are characterized by the different gain factors between the 8 and 13 TeV cross sections and by different signal acceptances for the experimental searches, possibly corresponding to different event selection categories. As can be seen from the numerical values in table 4, the geometric acceptances for the various production channels are quite similar to each other. The most important differences, of the order of ∼ 20%, are present for the CMS analysis, which explicitly presents the results in two categories: barrel-barrel (EBEB), which includes events with both photons in the central detector region, and barrel-endcap (EBEE), in which one photon is central while the second falls in the detector endcap. As we discussed before, the various production channels lead to slightly different rapidity distributions for the final photons, thus giving rise to different acceptances for the two CMS categories. Obviously this property cold be used to differentiate the production channels, although at present the experimental sensitivity is limited. We will discuss better this aspect in Appendix B.1.
Let us start the description of the numerical results by considering the quark and gluon production modes. The difference between the gg mode and the production through sea quarks (ss, cc, bb) is very small. All these channels have comparable gain factors (see table 3) and similar signal acceptances (see table 4). This is not unexpected, since the parton luminosities for these production channels are quite similar (see fig. 3). For this reason in our recast we only consider the gg channel, which provides also a good approximation of  Table 4: Acceptances for the scalar resonance case. The numerical results are derived for the following analyses: ATLAS 13 TeV [14], CMS 13 TeV (split into the two categories barrelbarrel (EBEB) and barrel-endcap (EBEE)) [15], ATLAS 8 TeV [18] and CMS 8 TeV [19].
The efficiencies for the gg case also apply to the ss, cc, bb, tt, since the differences among all these cases are 2%.
the sea quark ones. Significant differences, instead, are present with respect to the valence quark modes (uu and dd), mostly due to the gain factors that are much smaller than in the gg case. In addition to the acceptances we also included some reconstruction efficiency factors for the signal, which we take from the experimental papers. The numerical values are 70% for ATLAS 13 TeV, 81% and 77% for the EBEB and EBEE categories of the CMS 13 TeV analysis, 56% for ATLAS 8 TeV and 81% for CMS 8 TeV. Finally, in our numerical analyses we assume the resonance to have a small width, below the experimental resolution ∼ 7 GeV. The local significance of the diphoton excess is shown in the left panel of fig. 4 as a function of the R gg , R uu and R dd parameters. Since these tree parameters sum up to one, it is convenient to present the results in a "triangle" plot. One can see that the local p-value is sensitive almost exclusively to R gg , ranging from 4 σ in the case with purely gg-initiated production (R gg = 1) to 3.5 σ in the cases with R gg = 0. The dependence on the other two parameters is quite limited, since the gain and efficiencies for the uu and dd modes are similar. The best fit of the signal cross section is shown in the right panel of fig. 4 and ranges from 5 fb for the R gg = 1 case to 3 fb for R gg = 0.
In fig. 5 we show the combined goodness of fit (see Appendix B for more details). One can see that the compatibility of the various searches is never high. In the best case R gg = 1, the compatibility is only ∼ 9%, while it drops below 1% in the uu and dd-initiated modes. Analyzing the breakdown of the likelihood in each experimental search, one finds that the major source of tension is the ATLAS 13 TeV search, which favors a quite large signal cross sections ∼ 10 fb, to be compared with the much smaller ones ∼ 2 fb preferred by the other three searches. On the other hand, the two CMS searches and the ATLAS 8 TeV one show a very good compatibility ( 30%). Obviously the better global agreement found in the case of gg-initiated production is due to the larger gain factor between the 8 TeV and 13 TeV cross sections.
As a final case we consider the scenario in which the scalar resonance is produced exclusively through the γγ mode. 13 In this case there is no free parameter and the scenario is fully characterized by the gain factor r γ = 2.9 and by the efficiencies given in table 4. The efficiencies are quite similar to the ones for the dd initiated mode, thus we expect the overall features of this scenario to be comparable to the case R dd = 1. As we already discussed, the 13 For phenomenological analyses of this scenario see ref. [20].  relatively small gain factor increases the tension between the ALTAS 13 TeV results and the other searches, thus the photon-initiated mode is not favored by the present data. The goodness of fit is indeed 1% and also the statistical significance of the excess is relatively small, 3.5 σ. The best fit of the 13 TeV signal cross section is 3 fb.

Spin-2 resonances
Let us now move to the case of spin-2 states. As we did for the scalar resonances, we will adopt here a broad perspective and we will consider a generic new state without imposing any restriction on its production modes and on its decay distributions. When looking for a physics interpretation of these scenarios, it must be however kept in mind that resonances of spin J ≥ 2 have a typical interpretation as composite states. Therefore, the hypothesis that a new state of this kind is within the reach of the LHC requires an exotic strong dynamics not far above the TeV scale. Such a framework is considerably constrained by a variety of experimental tests, which limit the number of realistic benchmark scenarios. As for the scalars, the production modes can be encoded in the R in parameters defined as in eq. (19). In the present case, however, additional free parameters are needed to take into account the angular distribution of the decay products. As we explained in section 2, the decay distributions in the COM frame are a combination of a limited number of functional forms, which depend on the production channel (three forms for the gg and γγ mode and four for the quark-initiated channels). The total number of free parameters is thus significantly greater than in the scalar-resonance case. It is thus unpractical to keep all of them free in an analysis, but instead it is reasonable to consider a few benchmark scenarios. In the following we will describe some of them. In particular we will focus on single production modes, namely the gg initiated channel and the quark production modes. In addition we will also discuss a benchmark that parametrizes a very specific, but well motivated scenario,  the Randall-Sundrum graviton. 14

gg-initiated production
The most straightforward way to couple an exotic strong dynamics to the SM is via gauge interactions. This is typically realized whenever the constituents of the resonance are charged under the SM gauge symmetry. If the strong sector is charged under QCD, the leading production modes at hadron collider is expected to be the one involving gluons. Another interesting possibility is the case in which the resonance is produced from photons. Applying the results of section 2, it is straightforward to check that the COM angular distribution of the decay products is a combination of the functions D 2,0 (see table 2), we are left with just three possible functional forms. We can thus fully parametrize the differential cross section as 0,0 P 00 + D   as a function of three free quantities, P 00 , P 02 + P 20 and P 22 , which are normalized such that they sum up to unity.
As a representative example, we recast the experimental searches for a diphoton resonance in the scenario with a narrow spin-2 resonance produced exclusively from gg. The case of γγ production is similar, however, analogously to the scalar case, it is disfavored by the current data because of the small cross section gain between 8 and 13 TeV.
The geometric acceptances for the various experimental searches are listed in table 5 (see  table 6 for the acceptances in the γγ channel). In fig. 6 we show the signal significance and the best fit of the cross section for the gg mode as a function of the three free parameters, P 00 , P 02 + P 20 and P 00 . The goodness of the fit is instead shown in fig. 7. We find that the signal significance is around 4 σ and is slightly higher for a resonance decaying in the D 2,0 modes. The goodness of fit in the P 02 +P 20 = 1 corner is ∼ 12% and is significantly higher that in the other configurations, in particular for P 22 = 1 we find a compatibility around 4%. The best fit of the signal cross section varies from ∼ 4 fb in the configurations with P 02 + P 20 = 1 to ∼ 7 fb in the cases P 00 = 1 and P 22 = 1.

qq-initiated production
A spin-2 resonance can have sizable couplings to quarks if some of the latter mix significantly with fermionic composites of the exotic strong dynamics. We can thus envisage a scenario in which a spin-2 resonance is produced mainly through the qq channel. In this set-up the initial partons can have m = ±1 or m = 0. The latter spin, however, is only generated by interactions suppressed by a chirality flip, which thus are expected to give rise to smaller contributions than the |m| = ±1 channel. For this reason we will neglect the m = 0 case in what follows. In the m = ±1 channel, the decay distribution can be parametrized in terms of two quantities, P 10 and P 12 , so that 1,0 P 10 + D 1,2 P 12 .
In principle all quarks could couple to the new resonance. However in order to avoid tensions with existing bounds we assume the resonance has negligible flavor-violating couplings to the light quarks. Plausible scenarios may be constructed if the coupling is either family-universal or dominantly with the heavy quarks (in particular with the third generation). In the following we will thus consider two benchmark scenarios. In the first the spin-2 resonance couples dominantly to the bottom quark. In the second scenario it has a family-universal coupling with a single quark representation (as, for instance, the righthanded up-type quarks). We will assume that in both scenarios the unavoidable coupling to gluons can be neglected.
The geometric acceptances for the various quark production channels are listed in table 7. The signal significance and the best fit of the cross section for the bb production channel is shown in fig. 8 as a function of the P 10 = 1 − P 12 parameter. One can see that the significance is around 4 σ and the cross section best fit is 5.5 fb. This scenario provides a good compatibility among the experiments, at the level of 10−15%. The dependence on P 10 is relatively mild due to the limited statistical precision currently available. As can be seen   Table 8: Acceptances for an RS1 graviton into photon pairs. from table 7, the acceptances for the two different angular distributions differ significantly, thus they could allow to better differentiate the various scenarios when more data will be available.
The results for the scenario with universal couplings to the fermions are shown in fig. 9. In this case the cross section gain factor is mostly determined by the one of the valence quarks and is given by r univ = 2.9. Due to the relatively small gain factor the significance in this scenario is lower than in the bb channel, namely it is around ∼ 3.5 σ. The best fit for the signal cross section is ∼ 3 fb. In this scenario the compatibility among the experiments is rather poor, at the 1% level.

The RS graviton
We conclude the discussion of the spin-2 resonances by considering a well-known scenario that includes a new state of this kind, the Randall-Sundrum (RS1) model. In this case the spin-2 state is identified with the massive graviton of RS1, which is coupled to the stress tensor of the SM. The particular form of the coupling implies a peculiar relation between qq and gg production modes, namely σ qq = σ qq = 2 3 σ gg . It also completely fixes the angular distributions of the diphoton final state. In the notation introduced in the previous section one gets P gg 22 = P qq 12 = 1. On top of fixing the properties of the diphoton final state, the RS1 scenario also determines the relative importance of the other decay channels of the massive graviton. In particular it possesses large branching ratios into leptons, which suggests that diphoton searches are probably not competitive with di-leptons for this specific scenario. Taking into account the implications of the other final states, however, goes beyond the scope of this paper, thus we just concentrate on the diphoton channel. The geometric acceptance for the RS1 graviton are listed in table 8, while the gain factor between the 8 and 13 TeV cross section is mostly determined by the gg production mode and is equal to r RS = 4.1. We find that the available searches imply a statistical significance of 3.8 σ for a graviton with a mass 750 GeV, with a compatibility of the different searches at the 3% level. The best fit of the signal cross section is 5 fb.

Spin-3 resonances
As a last benchmark scenario we discuss the case of spin-3 resonances. From table 1, one sees that for a resonance of odd spin produced through gg or γγ only the channels |m| = 2 are allowed. Therefore the most general gg/γγ cross sections can be written in terms of a single parameter dσ in d cos θ =σ in D   Table 9: Acceptances for a spin-3 resonance produced in the gg channel.
Quark production, analogously to the even-spin case, can instead occur via both m = 0 and |m| = 1: Here, for simplicity, we will focus on the case of a spin-3 resonance produced in the gg channel. The acceptances for this scenario are listed in Table 9. From our recast of the experimental searches we find that the hypothesis of a resonance with a mass of 750 GeV has a statistical significance of 4.2 σ with a best fit of the signal cross section 5.6 fb. The compatibility of the various experimental results is 14%. The higher significance and better compatibility between the various searches comes from the fact that the decay distribution of the two photons (controlled by D 2,2 ) is quite central (similarly to D (2) 0,2 for the analogous spin-2 benchmark model). This implies a larger geometric acceptance for the ALTAS 13 TeV search and a slightly lower acceptance for the other searches. This difference mitigates the preference for higher signal strengths implied by the ATLAS 13 TeV data.

Conclusions and outlook
In this paper we provided a general characterization of the resonant diphoton production at hadron colliders. Our main result is the derivation of a new, simple phenomenological parametrization that can be used to describe resonances with arbitrary (integer) spin and CP parity, produced in any of the gg, qq and γγ partonic channels. By exploiting angular momentum conservation, the decay distributions of the resonance can be expressed as a combination of a small number of basis function that encode the angular distributions of the diphoton pair in the COM frame. The form of the basis functions is fully determined by the spin of the resonance. Their relative importance in the signal distributions, as well as the relative importance of the various production channels, are controlled by polarized resonance cross sections and decay branching ratios.
An important advantage of our parametrization is the fact that it does not depend on any assumption about the underlying theory describing the resonance. In particular it can be used even if the resonance dynamics can not be encoded into a local effective Lagrangian, which could be the case if it emerges from a strongly-coupled QCD-like dynamics. Our approach is thus completely model-independent and particularly suitable to describe in an unbiased way a possible signal observed in the diphoton channel. Although mainly aimed at characterizing a possible signal, our parametrization can also be used to express exclusions in the case of a measurement compatible with the background-only hypothesis.
As an example of the use of our results, we performed a simple recast of the ATLAS and CMS resonant diphoton searches, which recently reported an excess around an invariant mass M = 750 GeV. These recasts should not be interpreted as fully quantitative results, but rather as an illustration of the usage of our parametrization. For definiteness we focused on a few benchmark scenarios with resonances of spin J = 0, J = 2 and J = 3.
The J = 0 case is particularly simple, since the diphoton angular distribution is fixed to be completely flat in the COM frame. The properties of the resonance thus only depend on the relative importance of the various partonic production channels. Each channel is characterized by the gain ratio between the 8 and 13 TeV production cross section and by the acceptances in the various searches, which depend on the y distribution. We found that the gg, ss, cc and bb channels are quite similar and difficult to distinguish experimentally. The situation is instead different for the uu, dd and γγ channels, which have a significantly smaller gain ratio with respect to the gg mode. The present data show some degree of tension between the 8 TeV results and the 13 TeV ones, in particular the ATLAS analysis, which prefers large gain ratios. As a consequence the gg or heavy-quarks production modes are favored. In this cases a good signal significance, ∼ 4σ, is found with a 9% compatibility among the various searches. The compatibility is instead poor, around 1%, in the case of the light-quarks or γγ production modes.
Since they can lead to different non-trivial angular distributions for the diphoton pair, spin-2 resonances are characterized by a more varied phenomenology. Also in this scenario production channels, as the gg one, with large gain ratios are preferred. Moreover more central angular distributions are slightly favored since they lead to a higher acceptance, especially in the 13 TeV searches (see table 5). In the most favorable case, namely gg production with the D (2) 0,2 angular distribution, a signal significance of ∼ 4.2σ is found with a compatibility of 12% between the various searches. Another scenario that has a good compatibility with the data is the case of bb-initiated production, which can lead to an overall compatibility of 15%. Another spin-2 benchmark we considered is the case of a Randall-Sundrum massive graviton. In this scenario the production mode is dominantly gg and the angular distribution is described by the function D 2,2 , which leads to a more forward diphoton distribution. This property implies a not so good compatibility with the data, at the 3% level.
As a final scenario we considered a spin-3 resonance produced in the gg channel. This set-up is particularly simple since it is characterized by a single angular function, D 2,2 . In this case we find a good significance ∼ 4.2σ and a good compatibility among the various searches, ∼ 14%.
Besides providing the general framework within which benchmark scenarios can be defined, the phenomenological analysis presented in section 2 allows us to draw interesting conclusions concerning which properties of the resonance (once it is discovered) could be extracted from a careful experimental study of the resonant diphoton signal. Namely, we saw that the resonance spin and production mode could be established, barring peculiar degeneracies which we have identified, from the combined measurement of the cos θ and y distributions. Within a given hypothesis for the resonance spin and production mode, the cos θ distribution also gives us information about the resonance CP-parity. Indeed nonvanishing A ±∓ amplitudes (recall that the a's in table 1 are CP-even while the a's are CP-odd), which we could detect through the presence of a D 1,S or D 2,S component in the angular distribution, would imply either that the resonance is CP-even or that CP is badly broken by the resonance couplings. If instead A ±∓ were to vanish, we would not be able to distinguish a CP-odd R from a CP-even resonance with accidentally vanishing a 1,−1,2 . The only way to achieve this would be to measure a 0 and a 0 separately, but this is impossible since only a combination of the two enters, through eq. (11), in the differential cross section. This problem is particularly severe for J = 0, where A ±∓ = 0 by spin conservation and thus the resonance CP-parity cannot be measured.
A possible way out is to study, as pointed out in refs. [2,22] for the J = 0 case, the structure of the forward initial state radiation (ISR) that unavoidably accompanies the hard resonance production process. Consider the emission of two forward ISR jets 16 emitted in the forward and backward direction, respectively, and denote by ϕ 1 and ϕ 2 their azimuthal angles. For p ⊥ (j 1,2 ) M , the Feynman amplitude for the complete 2 → 4 process takes the form [23] A(in → j 1 j 2 γγ) ∝ where φ is the azimuthal angle of the hard scattering plane, i.e. the one of the diphoton pair appearing in eq. (3). The e iλϕ factors from the parton splittings are dictated by momentum conservation, as discussed in ref. [23] for the case of effective massive vector bosons splittings. The g λ 1,2 's are given functions, specific of the ISR splitting process at hand, of the incoming partons momentum fractions x 1,2 . The above formula illustrates that by studying the kinematical distributions of the ISR jets one can get more information about the polarized resonant production amplitudes than that obtainable from the 2 → 2 process. Taking for simplicity the soft limit, in which the most singular g-functions (i.e., those for gg, qg, and qq splittings) become independent of λ, one easily obtains an approximate formula for the complete 2 → 4 process cross section. Such a cross section, differential in the azimuthal angular difference between the two jets, ϕ 12 = ϕ 1 − ϕ 2 , and integrated over all other variables, reads for, respectively, gg/γγ and qq hard production. We defined δ = arg(A ++ /A −− ). Because according to table 1 a CP-even (CP-odd) resonance has δ = 0 (δ = π), we see that measuring the ϕ 12 distribution one would be able to infer the resonance CP-parity, even for J = 0. The distribution can also tell us if the a's in table 1 are complex, which would mean that the resonance interactions are mediated by loops of light particles as we discussed around eq. (5). Indeed, it might allow to extract the ration |A ++ |/|A −− |, which is necessarily equal to one if the a's and a's are real. However |A ++ |/|A −− | = 1 is also ensured by the CP symmetry, therefore observing |A ++ |/|A −− | = 1 would also mean that CP is broken. Another process which is worth considering, because of its larger rate, is the emission of a single detectable forward jet, with azimuthal angle ϕ j . In this case one must study the doubly differential distribution in cos θ and in ϕ = ϕ j − φ, i.e. the angle between the jet and the diphoton plane. The angular dependence, focusing once again on the soft/collinear limit and assuming for simplicity a heavy mediator (real a's), is controlled by for gg/γγ hard production, and for qq (a similar result holds for qq). Here the d's are the Wigner matrices for a generic spin J and BR S is the polarized branching ratio of eq. (12). Differently from the 2-jets emission previously discussed, studying the single ISR jet distribution does not furnish conclusive information about the resonance CP-parity at J = 0 because the dependence on ϕ disappears in the scalar case. Still, the measurement of this process gives access to different parameter combinations which do not appear in the fully inclusive 2 → 2 reaction and thus it is nevertheless worth studying. A detailed analysis of the ISR radiation pattern, and its potential implications for the experimental characterization of the resonance properties, is left for future work.

A On-shell amplitudes
In this appendix we will derive the effective couplings that parametrize the on-shell dynamics of a spin-0 or spin-2 resonance decaying into a photon pair. We will first compute the on-shell amplitudes for the production of R, from which the A in λ 1 λ 2 immediately follow. The analogous amplitudes for R → γγ can be straightforwardly obtained from them. In section A.3 we will then present an effective Lagrangian that may be employed to implement the relevant processes into a Montecarlo generator.
The amplitudes for in → R depend only on a few basic quantities. First of all they are a function of the 4-momenta of the initial partons, which we denote by p µ 1 and p µ 2 . For later convenience, we introduce the notation p µ = p µ 1 + p µ 2 for the resonance momentum and q µ = p µ 1 − p µ 2 for the other independent combination of the initial momenta. The only non-trivial Lorentz scalar is given by the resonance mass, p 2 = −q 2 = M 2 . The amplitudes also depend on the polarization vectors µ 1,2 (for gg/γγ production) and the spinors u 1 and v 2 of the SM quarks (for qq production). In the case of a spin-2 resonance, an additional tensor t µν is present, that describes the polarization of R.
Since the external states are on-shell, the equations of motion may be used to simplify the expressions. Specifically, the polarization tensor t µν of a spin-2 resonance is required to satisfy the same conditions as an on-shell spin-2 field, i.e. to be transverse and symmetrictraceless. This implies that its contraction with p µ vanishes and the only non-trivial terms can be obtained by contractions with q µ and µ 1,2 . Similarly, the equations of motion for the SM quarks can be used to remove factors of γ µ p µ i from the amplitudes relevant to qq production. Furthermore, because of the transversality of the gauge bosons µ i p i µ = 0 may be used to simplify the expressions for γγ, gg production and γγ decay. Finally, Lorentz invariance constrains the form of the amplitudes for gg, γγ → R (and analogously R → γγ) in a non-trivial way, by forcing them to be invariant under a shift µ i → µ i + p µ i , where p µ i is on-shell [24]. Amplitudes consistent with Lorentz invariance therefore automatically satisfy the on-shell Ward identities.
Interestingly, we find that the amplitudes for gg, γγ → R derived following the above recipe can be unambiguously uplifted to expressions valid for off-shell gauge bosons (and not necessarily transversely polarized) and respecting the off-shell Ward identities. This allowed us to perform a sum over gauge boson helicities in the familiar fashion µ ( ν ) * → −g µν and check that the squared amplitudes thus obtained agree with the ones derived from the helicity amplitudes and the general formalism of section 2. For completeness we will present our results in this off-shell form.
In the following we will specialize the discussion to the J = 0 and J = 2 cases, although resonances with higher spin can be treated analogously. Our results agree with ref. [8] up to phase conventions.

A.1 Spin-0 resonance
As a first case we consider the production amplitude for a spin-0 resonance. In the gg or γγ channels we find For in = qq instead we obtain In the above equations, the CP-even (CP-odd) coefficients a 0 , b 0 are all dimensionless quantities (in general complex). The helicity amplitudes A λ 1 λ 2 can be straightforwardly derived from eq.s (25) and (26). For the polarization tensors (including both J = 1, 2) we use the conventions of [25], whereas the spinors are taken from [26]. The result reads in agreement with table 1.

A.2 Spin-2 resonance
We can now discuss the scenario with a spin-2 resonance. In this case there are two differences with respect to the spin-0 state. First, the amplitude will depend on the polarization tensor of the resonance t µν . Second, the amplitudes for gg/γγ → R, which according to the rules described above may contain terms that are anti-symmetric in the exchange of the incoming state particles, must be symmetrized since the two gauge bosons are indistinguishable. Starting with gg, γγ → R, we obtain To arrive at these expression we used the Schouten identity to eliminate all CP-odd structures in which t µν is contracted with µ i or with the Levi-Civita tensor. These either vanish identically or are equivalent to a renormalization of the vertices in eq. (28). In the case in = qq we find From eq.s (28) and (29) one can derive the corresponding helicity amplitudes again in agreement with table 1.

A.3 On-shell Lagrangian
Finally we present two effective Lagrangians that may be used to simulate spin-0 and spin-2 diphoton resonances through Montecarlo generators. The various terms appearing in the production amplitudes in eq.s (25), (26) and (28), (29) may be thought of as effectively arising from the following set of effective operators + R −a q 0 qq + i a q 0 qγ 5 q , in the case of a spin-0 resonance, and in the case of a spin-2 state. In our notation F µν is the field strength of either photons or gluons and F µν = 1 2 µναβ F αβ . We warn the reader that eq.s (31) and (32), where in general a, a are complex , are not meant to describe the off-shell dynamics of these models. They just represent a practical parametrization of the on-shell couplings relevant for resonant production. Even if R ultimately arises from a theory consistent with Lorentz and gauge symmetries that does not admit a Lagrangian formulation, the effective description (31) (32) can be used to parametrize the diphoton resonant production and decay.

B Statistical treatment
In this appendix we briefly describe the statistical procedure we used to obtain the numerical results for the benchmark scenarios presented in section 3. For our recast we followed a simple strategy to reconstruct the profile likelihood ratio of the various searches (q(µ), where µ is the signal strength parameter) by exploiting the available data, namely the p-value of the background-only hypothesis and the exclusion limits on the signal cross section.
The p-value is directly connected to the value of the profile likelihood ratio for a vanishing signal hypothesis q(µ = 0). In the asymptotic limit, q follows a half-χ 2 distribution [27], and the background-only p-value corresponds to the cumulative distribution starting at q(µ = 0). 17 From the exclusion limits, instead, one can reconstruct the value of the cross signal strength µ for which the cumulative distribution is equal to the exclusion threshold. These elements are enough to reconstruct the profile likelihood if we assume that a Gaussian approximation is valid. In this case the profile likelihood is just a quadratic polynomial in µ, that is always positive and vanishes in a single point, i.e. for µ equal to the signal strength best fitμ. 18 Notice that with our procedure we are only able to extract the global likelihood ratio for each experiment, but we have no access to the likelihood for the single event categories used in the experimental analysis. For instance in the CMS 13 TeV analysis two categories are considered which could give some information on the angular distribution of the diphoton events. Due to the limited statistics, however, this information is not extremely significant and our approximate results are reliable. We will discuss this point quantitatively in subsection B.1. The above procedure can be straightforwardly applied to the CMS 8 and 13 TeV analyses and to the ATLAS 13 TeV one, which provide the p-value and the exclusion limits as a function of the diphoton system invariant mass for the case of a narrow-width resonance. In the recast it is important to take into account the fact that the CMS results are provided for a scenario with a RS graviton, while ATLAS consider the case of a scalar resonance. This implies different acceptances and reconstruction efficiencies as we discussed in section 3. For the ATLAS 8 TeV search, on the other hand, only the exclusion limits are publicly available. In this case we reconstructed the profile likelihood by estimating the width of the 1 σ and 2 σ bands from the expected exclusion limits. Since in this search the data show only a very mild (below 1 σ) excess, our estimate is expected to be fairly accurate. We also checked that this procedure is correct by using the CMS 8 TeV data, in which case we find that the reconstructed likelihood is very close to the one obtained with the other method we used.
Once the likelihood ratios for the various searches are reconstructed, it is straightforward to use them to extract the best fit of the signal strengthμ and the combined signal sig- 17 Notice that the ATLAS collaboration in the analysis of the 13 TeV data used a slightly different procedure, the uncapped p-value. This definition, however, coincide with the usual one if the best fit of the cross section is for µ > 0. This is always the case if there is an excess in the data, as it happens in the available ATLAS and CMS 8 and 13 TeV searches for mγγ ∼ 750 GeV. 18 This procedure is correct in the case in which an excess is present in the data, in which case necessarilŷ µ > 0. In the case of a deficit of events the signal strength best fit would be negativeμ < 0, but the profile likelihood is defined in such a way to vanish for µ = 0. In this case the knowledge of the background-only p-value and of the exclusion limit is not enough to fully reconstruct the likelihood ratio.  nificance, i.e. the p-value of the background-only hypothesis. Another interesting quantity that can be computed is the compatibility among the various searches, also known as the "goodness" of the fit [28]. To extract this quantity one compares the likelihood for the best fit of the cross section with the one obtained by assuming independent signal strengths for each experimental search. The resulting likelihood follows a χ 2 distribution with a number of degrees of freedom equal to the number of experiments minus one.
As an application of our recast procedure we show in fig. 10 the statistical significance of the signal for the scenario with a narrow scalar resonance produced in the uu, gg and bb channels. In the plots the p-values for the single searches are shown as a function of the resonance mass, together with the result for the combination of the 13 TeV searches only and the full combination of the 8 and 13 TeV data. One can see that the significance of the full combination for M 750 GeV is quite close to the one of the 13 TeV only searches if the resonance is produced in channels with a large cross section gain between 8 and 13 TeV, namely the gg and bb modes. This shows that in these scenarios the agreement between the 8 TeV and the 13 TeV data is reasonably good. On the contrary, in the uu case, the p-value for the full combination is significantly smaller than the one for the 13 TeV searches only, implying a sizable degree of tension among the experimental searches. These results confirm what we found in section 3.
Finally in figs. 11 and 12, we provide the fit of the signal for the scenario with narrow B.1 Impact of the CMS 13 TeV categories As we saw in the section 3, the impact of the angular distribution on the various searches can be significant, even if we only consider the total number of events without explicitly looking at the distributions. The reason for this dependence is the fact that relatively hard cuts are imposed on the signal, with the aim of selecting events in which the final-state photons are central. As a consequence angular distributions that enhance the signal in the central region of the detector have larger acceptances than the ones that give rise to a more forward signal. More details on the angular distribution can in principle be obtained by looking at the different signal categories used in the experimental analyses. In particular the CMS 13 TeV study splits the events in two categories: the EBEB in which both photons are in the barrel of the detector (|η| < 1.44) and EBEE in which one photon is in the barrel while the second is in the endcap (|η| ∈ [1.57, 2.5]). 20 The two categories allow to get a rough information on the angular distribution, thus improving the discrimination power for the different benchmark models. Unfortunately, the procedure we described in the previous section to recast the experimental analyses did not allow us to take into account separately the different categories. We now want to estimate how drastic this simplification is and how much a full analysis could help in discriminating the angular distribution of the diphoton signal.
For this purpose, here we implement a simple recast of the CMS 13 TeV search. 21 We reconstruct the likelihoods associated to the two event categories by using the distributions of events provided in fig. 3 of ref. [15]. We assume that in each bin the events follow a Poisson distribution and we construct the total likelihood by multiplying the likelihoods for each bin. We model the background using the functional form given in the experimental paper f 0 (m γγ ) = N e −p 1 mγγ m −p 2 γγ , where p 1,2 are free parameters and N is the overall normalization which we fit together with the other parameters. For simplicity we only focus on the narrow resonance scenario, and we model the signal by a Gaussian distribution with a half width equal to the experimental resolution (10 GeV for the EBEB category and 16 GeV for the EBEE category). The signal and the background are fitted simultaneously for each signal strength hypothesis. The test statistics we use is based on the profile likelihood ratio and the background-only p-value is computed by assuming that the distribution is asymptotically equal to a half-χ 2 distribution with one degree of freedom. The result of our recast is shown in Fig. 13, where we plot the significance of the signal as a function of the acceptance in the EBEB category Acc EBEB under the assumption that the total acceptance is equal to 1. One can see that the statistical significance of the signal has a non-negligible dependence on the angular distribution. In particular the signal is mostly present in the EBEB category, so that models with a more central signal distribution are The results are shown as a function of the relative signal acceptance in the EBEB category (Acc EBEB ), while the acceptance in the EBEE category is Acc EBEE = 1 − Acc EBEB . The unshaded region denotes values of Acc EBEB that can be obtained in the various benchmark models we discussed in the main text.
preferred. The present experimental sensitivity, however, is not very large, so that the impact on the fit in the benchmark models we considered is mild. This justifies our approximation of combining the two CMS 13 TeV categories. 22