#### Abstract

LHC run-II has a great potential to search for new resonances in the diphoton channel. Latest 13 TeV data already put stringent limits on the cross sections in the diphoton channel assuming the resonance is produced through the gluon-gluon fusion. Many beyond the Standard Model (SM) theories predict TeV-scale scalars, which copiously decay to diphotons. Apart from the gluon-gluon fusion production, these scalars can also be dominantly produced in other ways too at the LHC, namely, through the quark-quark fusion or the gauge boson fusions like the photon-photon, photon-, , or fusions. In this paper we use an effective field theory approach where a heavy scalar can be produced in various ways and recast the latest ATLAS diphoton resonance search to put model-independent limits on its mass and effective couplings to the SM particles. If a new scalar is discovered at the LHC, it would be very important to identify its production mechanism in order to probe the nature of the underlying theory. We show that combining various kinematic variables in a multivariate analysis can be very powerful to distinguish different production mechanisms from one another.

#### 1. Introduction

There are numerous theoretical motivations to expect that the Standard Model (SM) is not the complete story and the scalar spectrum of a larger theory may be richer than to possess only one neutral scalar, the Higgs boson. From this expectation searches for new scalars are continuously being carried out at the Large Hadron Collider (LHC) in various channels. Along these directions, no confirmed hint has been found so far in any of these searches. Nevertheless, some anomalies in some of these searches have drawn significant attention in the high energy physics community, recently. Among them, the most famous one is the 750 GeV diphoton excess [1, 2] which created a lot of excitements in the community. Before the excess went away with more data, numerous attempts have been made to explain the excess (see [3] for a review and a long list of references). Another important excess was the diboson excess around 2 TeV resonance mass [4–6] which also later turned out to be a statistical fluctuation.

Searches for heavy scalars at the LHC are generally being carried out in the diphoton, diboson, or dijet resonance searches. The diphoton channel, among them, is particularly important as this channel provides a comparatively cleaner background. Higgs boson was first discovered in the diphoton channel at the LHC [7, 8]. TeV-scale scalars decaying into a diphoton system is one of the key predictions of many beyond the Standard Model (BSM) theories. Various possibilities have been extensively explored in the context of the 750 GeV diphoton excess (see the reference list of [3] for various models that predict TeV-scale diphoton resonances). To test these predictions, LHC run-II provides us with a great opportunity to observe a diphoton resonance of mass up to a few TeV. In this paper we particularly focus on the diphoton final state for these reasons.

If a particle decays to diphotons, it must either be spin-0 or spin-2 in nature as the spin-1 particles decay to on-shell diphotons is forbidden by the Landau-Yang theorem [9, 10]. A spin-2 particle or graviton couples universally to all matter fields through energy-momentum tensor. Various extra dimensional models like the ADD model [11] or the RS model [12] predict the existence of graviton. If a resonance in the diphoton system mediated by graviton is observed, one would expect resonances at the same mass in other possible channels also. Therefore, simultaneous studies in various channels might be more illuminating for the spin-2 particle. The current limits on the graviton mass are already quite high, around ~2-3 TeV [13, 14]. On the other hand, scalars of mass ~1 TeV decaying into diphotons, which is a typical signature of many models, are still allowed by the LHC data. These scalars can be produced at the LHC in various ways, namely. through the , , , , , or fusions. In this paper we consider a model-independent effective field theory (EFT) approach where the scalar can be produced and decayed (two-body) in different possible ways as mentioned. But we only concentrate on the diphoton decay mode in this paper as stated earlier.

First, we derive the available parameter space for a scalar (produced in different ways) decays to diphotons using our EFT approach. These limits will be grossly model-independent and can be used to set limits on other models wherever applicable. If a scalar resonance will actually be seen in future, the most obvious question that will arise is how the scalar is produced? A most common way to decipher the production mechanism of a heavy scalar is to look at various kinematic distributions especially various jet observables, which are important in this regard. This has been investigated to some extent in the literature in the context of the 750 GeV diphoton excess [15–19]. In this paper we revisit some of the jet observables and show their effectiveness in distinguishing different production modes. We, then, use a multivariate analysis (MVA) by combining many kinematic variables to distinguish different production modes more efficiently.

This paper is organized as follows: in Section 2, we employ an effective Lagrangian for the scalar ; in Section 3, we discuss about the decays and various production modes of at the LHC and derive exclusion limits on the mass and couplings from the latest diphoton resonance search data. In the same section, we discuss how two different production modes of can be distinguished using a MVA analysis. Finally, we conclude in Section 4.

#### 2. Effective Lagrangian

We consider an EFT approach where a heavy scalar interacts with the SM gauge bosons through the dimension-5 operators and with the SM quarks through the dimension-4 operators. Assuming is a CP-even real scalar, we employ the following effective Lagrangian:where the field-strength tensors corresponding to gluon (), photon (), , and bosons are , , , and , respectively, and their generic form is where . All the dimension-5 operators are suppressed by the new physics scale . In general, could be different for different operators, but we assume they are the same for all the operators. Note that only the is a dimension-4 operator and we introduce the electroweak symmetry breaking scale GeV in association with to bring the scale in the interaction. The motivation behind this is that if operators are effectively originated from new physics then effective couplings are expected to contain the imprint of the scale (possibly in the form or some power of this ratio). This way of parameterizing the couplings also enables us to present exclusion limits on all the couplings in the form. Here, we use the notation to denote a generic dimensionless coupling associated with the vertex. The scalar can, in general, couple differently with the different SM quarks. For simplicity, in this analysis, we assume a single coupling , the same for all the SM quarks. Note that, in all interactions with the gauge bosons, the normalization factor is so chosen such that the corresponding Feynman rule takes the form where and are the 4 momenta of two gauge bosons and , respectively, directed towards the vertex. The Feynman rule for the interaction is .

In general, the new scalar can mix with the 125 GeV scalar () with a mixing angle . This leads to the scaling of all the couplings of by a factor . Although this would not change the branching ratios (BRs) of , it would change the production cross section of by a factor . Since all the measured signal strengths are pretty close to unity, this will make close to one. That is why, in this paper, we have neglected any mixing between and for simplicity.

#### 3. Phenomenology

In addition to the SM Lagrangian, we implement the effective Lagrangian of shown in (1) in FEYNRULES [20] to generate the Universal FeynRules Output [21] model files for the MADGRAPH [22] event generator. We use the MMHT14LO [23] parton distribution functions (PDFs) for event generation. This PDF set includes the photon PDF which has been computed following the approach described in [16, 24]. We use the factorization scale and the renormalization scale at in our analysis. Generated events are further showered and hadronized including multiple parton interactions by using PYTHIA8 [25]. We perform detector simulation using DELPHES [26] which uses FASTJET [27] for jet clustering. Jets are clustered using the anti- algorithm [28] with . We analyze the reconstructed objects by implementing ATLAS selection cuts [29], which we summarize in Section 3.3. For MVA, we use the adaptive Boosted Decision Tree (BDT) algorithm in the TMVA [30] framework.

##### 3.1. Decays of

From the Lagrangian in (1), we have the following two-body decay modes of , namely, , where . The partial widths for these decay modes are given by the following expressions: where denotes the electroweak gauge bosons and . There could be subdominant three-body decays of possibly mediated through an off-shell gauge boson. If the intermediate gauge boson is massless, in case of gluons or photons, the three-body BRs are nonnegligible especially when is large [31]. In this analysis, we consider the two-body and three-body decays of to obtain the total width where the three-body decay widths are computed numerically using MADGRAPH. Partial widths of three-body decay modes where an off-shell gauge boson goes to pair grow very rapidly with increasing scalar mass. This is due to the contribution coming from the longitudinal polarization of bosons. Therefore, in high mass region, BR for reduces substantially.

##### 3.2. Production of at the LHC

When all in (1) are nonzero, the scalar can be produced from the , , , , , and fusions at the LHC. In Figure 1, we show the partonic cross sections of different production modes of at the 13 TeV LHC for (taking one at a time) and TeV. In case of the production of through the or fusions, initial and come from the quark splitting. Therefore, is produced in association with at least two jets for this case. Similarly, for the initiated production, is produced in association with at least one jet. Partonic cross sections are computed by applying the following generation level cuts on the jets () and photons () wherever applicable Here, transverse momentum, pseudorapidity, and separation in the - plane are denoted by , and , respectively. These basic cuts are used to avoid any soft divergence present at the event generation level and stricter selection cuts are applied at the level of reconstructed event analysis after detector simulation. Note that all cross sections scale as , and therefore, we present them by choosing and TeV such that one can translate it easily for other values.

For all the six types of production of , we generate parton level events with up to two jets in the final state. These events are passed to PYTHIA8 [25] for showering and hadronization. This process may introduce double counting between the matrix element partons and the parton showers. To generate inclusive signal events by avoiding any double counting, we use the MLM matching [32] technique to match the matrix element partons with the parton shower. Inclusive signal events including up to two jets for the , , and fusions are generated by combining the following processes: where we set the matching scale GeV. The curved connections above two photons signify that they come from the decay of . To determine the appropriate for these production processes, we have done three important checks, namely, smooth transition in the differential jet-rate distributions between events with and jets; matched cross sections are within ~10% of the zero jet contribution and also do not vary much with the variation once we have chosen it properly. For the , , or fusion productions, the initial or come from the quark splitting and we have additional jets at the Born level process. Therefore, the and fusion events are generated only at the level and no matching is required for these cases. But for the fusion, we do use matching by combining the processes and with GeV. The dominant SM background (about 90% of the total) comes from the process. We generate this background by matching up to 2 jets with GeV.

##### 3.3. Exclusion from the LHC Data

Diphoton resonance searches at the LHC using run-I and run-II data set strong upper limits (ULs) on of a spin-0 or spin-2 resonances [29, 33]. It should be noted that these searches are generally optimized for an -channel resonance production through the fusion followed by its decay to two photons. If the resonance is not produced from the fusion, the selection cut efficiencies can vary depending on the different production mechanisms of the resonance. For a particular production mechanism, it can also vary significantly on the number of selected photons and jets. Therefore, in order to derive exclusion limits on the model parameters by recasting the limits on from an experiment, one has to properly take care of the selection cut efficiencies. This can be done properly by using the following relation [31]: where is the UL on the number of signal events, which can be written as the product of the signal cross section (produced through a particular mechanism used in the analysis), the corresponding signal cut efficiency , and the luminosity . When different types of production mechanisms contribute to any experiment, can be expressed by the sum . Here, runs over all the contributing production mechanisms.

To see the change in efficiency for the different production mechanisms and also for the different resonance masses, we roughly employ the following event selection cuts used by the ATLAS collaboration for their spin-0 diphoton resonance search as listed below [29].(1)Transverse energy of the two selected photons satisfies GeV and GeV and transverse momenta of selected jets satisfy GeV for and GeV for . Here, and denote the highest- and second highest- photons, respectively.(2)Pseudorapidity of the selected photons satisfies excluding the barrel-endcap region and jets .(3)Separation in the - plane between the two selected photons or any photon-jet or jet-jet pair satisfies .(4)Invariant mass of the two selected photons satisfies and .

In addition to the above set of cuts, we also apply default photon isolation cuts given in DELPHES for the ATLAS detector. Jets with high- mainly come from the vector boson fusion topologies. We use a threshold GeV for for better sensitivity [34]. In Figure 2, we show cut efficiencies for the cuts listed above for different production modes as functions of . The cut efficiency for ATLAS for their spin-0 resonance produced through the fusion is roughly about 62% [29] and we find very close agreement (around 60%) using our analysis codes. After validating our codes, we compute cut efficiencies for the other production modes for the selection cuts mentioned above and find that they do not vary much, only up to ~15% for different production modes. It is pointed out in the ATLAS paper [29] that the cut efficiencies for different production modes would not differ much for their signal criteria (fiducial region). As expected, in the high mass region, TeV, cut efficiencies become insensitive to the mass.

In our EFT approach, there are six free couplings that affect the production of . But taking all nonzero at the same time will make the analysis very complicated. Therefore, for simplicity, we choose only one as nonzero at a time, in addition to nonzero , and show the two-dimensional (2D) exclusion regions (colored) in the - plane for four benchmark masses, TeV (presented in Figure 3). Only in Figure 3(a), we show the exclusion regions (colored) in the - plane assuming all are zero except . To derive these limits, we recast 95% confidence level (CL) UL on the for the spin-0 resonance search by the ATLAS collaboration at the 13 TeV with [29]. This analysis is done for the resonance width MeV. If the width of a particle is very small compared to its mass, one can safely use the narrow width approximation (NWA). In all our results we use the NWA ignoring any interference effect between the signal and the background.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**The bumpy nature in the exclusion limit on in Figure 3(a) is due to nonsmooth nature of the observed UL on . The highest value of that is excluded is ~0.05 around TeV. By choosing a value for , one can translate this limit to . For instance, choosing for TeV, one finds that TeV is excluded. Basic shape of the exclusion regions in Figures 3(b) and 3(c) is different from the ones in Figures 3(d), 3(e), and 3(f). This is because the and the fusion productions dominate over the fusion mode for as seen in Figure 1. On the other hand, cross sections for the , , and fusion modes are smaller than the mode for unity . One should also notice that exclusion regions in Figures 3(d) and 3(f) are slightly asymmetric around axis. This is due to the interference effect between the and the or production modes. On the other hand, no interference is possible between the and the , , or fusion modes. Note that exclusion limits become insensitive to as we go to higher values. This is because the production cross section varies as and BR to diphoton varies as for region. This makes for large .

To include higher order effects, we use a constant next-to-leading order (NLO) -factor of 2 for the fusion [35]. The NLO corrections to a heavy scalar produced from the fusion are computed in [36] where it is found that the NLO -factor is close to 1 for heavier masses. If the scalar is produced from the light quark fusions, one might expect slightly bigger -factor. Here, we assume it to be 1 since it is not available in the literature. For the , , , and we assume it to be 1.3 [37]. The actual values of the -factors for different channels can be slightly different from the constant values we have used but they have very little effect on the exclusion limits.

In general, when and any one are nonzero, the production cross section can be expressed as where is the interference contribution and () in the* r.h.s.* is the production cross sections through the () fusion (see Figure 1). These cross sections as functions of mass have been computed numerically by interpolating cross sections points in the mass range 0.5–2.5 TeV. When more than one are nonzero, the combined selection cut efficiency, in general, depends on and . Including cut efficiencies in (7) (omitting the functional dependence on ), we get where , are the cut efficiencies for the pure and pure fusion production modes, respectively, and they are functions of only, whereas the combined efficiency and associated with the interference term are functions of , , and . We have seen that is mostly sensitive to but not to the couplings. Therefore, for simplicity we use for region and for region. Branching fraction in the channel can be expressed as where ’s are known analytically from (3). Finally, we derive exclusion regions in Figure 3 by using (8) and (9) in (6).

In Figure 4, we show combined cut efficiency for the coupling assumption (all other are zero). This combined cut efficiency should lie between the two individual efficiencies and according to its definition. One might also be interested to see the behavior of , and for different coupling assumptions. In Figure 5, we show these three quantities in the 2D plane for the two cases: and (other are set to zero).

**(a)**

**(b)**

##### 3.4. Distinguishing Different Production Modes

A common way to distinguish different production modes of a heavy scalar is to scrutinize various kinematic distributions especially the jet activities associated with the scalar. It was pointed out in [15–19] that the jet multiplicity () distribution could be very important in this regard. In Figure 6, we show the normalized distributions for various production modes of the scalar and compare them with the SM prediction. These distributions are obtained assuming TeV at the 13 TeV LHC with 50 integrated luminosity with the diphoton invariant mass () satisfying GeV, in addition to the set of cuts defined earlier. Our jet selection cuts are GeV for and GeV for . The dominant background contribution of about 90% comes from the SM process and, in this analysis, we only consider this as the background which we estimate from our simulation. The error bars associated with the background represent the statistical uncertainly only. In reality, various components of systematic uncertainties like the jet energy scale, jet energy resolution, and uncertainty in the luminosity must also be considered to obtain the total uncertainty [33]. But the systematic uncertainty becomes small compared to the statistical one when background distributions are obtained from data.

**(a)**

**(b)**

It is visibly clear that the different production modes display very different jet multiplicity distributions. The distributions for the and the fusion modes peak at 0 jets but the peak for the mode is not as sharp as the mode. Cross section for the 0-jet bin for the mode is roughly about 60% of the total cross section. On the other hand, it is about 45% for the fusion case. The SM background distribution also peaks at 0 jets but contains only 30% of the total cross section. The fusion shows a peak at 1 and 2 jets whereas the vector boson fusion production through the and fusions shows peak at 2 and 3 jets. The fusion mode, on the other hand, shows a peak at 1 jet.

Different nature of the distribution can be captured by the average jet multiplicity associated with the diphoton resonance. We compute the average jet multiplicity of different production modes and the background and report these numbers in Table 1. It is expected that if the scalar is produced through the fusion, then the average jet multiplicity is lower compared to the or the fusions. This is because, in case of the fusion, a hard jet in the final state can originate from the splitting. However, this is suppressed compared to the leading order (LO) process with zero jet by the small probability of branching and also by the smallness of . On the other hand, colored particles in the initial state, that is, in case of the or the fusions, lead to higher jet multiplicity. The average jet multiplicity is greater than two for the vector boson fusion modes because one would anticipate to get at least two hard jets in most of the events since two initial ’s come from the branching. For the initial state, one expects at least one hard jet from the splitting.

##### 3.5. Multivariate Analysis

In the previous subsection, we show as a demonstration that the jet multiplicity distributions of two different production modes (and also for the background) can be quite different. Apart from the distribution, there are other kinematic variables which also show some differences in their shapes for different production modes. For example, in [15], the authors showed that various distributions like the scalar sum of transverse energy , pseudorapidity () of the selected photons and jets, central rapidity gap () between the jets, and the scalar show some visible differences for the and the production modes.

A cut-based analysis which employs a set of rectangular cuts may not perform well to decipher the underlying production mechanism of the scalar. In order to effectively distinguish two different production modes, one can use various kinematic variables that show some (small) differences in their shapes simultaneously in a MVA whose output might show large differences in their shapes. If appropriate variables are chosen, a MVA is expected to perform better than a cut-based analysis. Generally, MVA techniques are used to separate signal from background. Here, we use a MVA technique (BDT) to distinguish two different production mechanisms more efficiently than a simple cut-based analysis. In particular, we use the adaptive BDT algorithm in the TMVA framework. We train the algorithm by tuning various parameters like the number of trees, minimum size of the node, and so forth for proper training of different production modes. Optimal values of these parameters are not fixed and they can differ for each analysis.

For MVA, we select events with at least one jet and construct twelve simple kinematic variables as shown in Table 2. This includes , , between and leading jet, and of two selected photons, and the leading jet and the separation in the - plane between the photons and the leading jet. These twelve variables are finalized from a bigger set of variables by looking at their discriminatory power and less correlation. In particular, the variables we use are not correlated more than ~40% for signal. But these correlations might be different for the background. Next to each variable in Table 2, we show their relative importance in the BDT response and these numbers are obtained from TMVA using the and the production modes. Relative importance is a fraction (where all importance sums up to unity) which is used to identify the ranking of the variables in MVA. In other words, greater relative importance of a variable signifies that the variable is a better discriminator. For actual definition of relative importance, interested readers may look into the TMVA manual. From Table 2, we see that is the best discriminator to differentiate the and the production modes. Other variables like , , , , and also act as good discriminators. Here, our main aim is to distinguish different production mechanisms of the scalar using a suitable MVA. Before arriving to this step, one might be interested to see comparisons of various kinematic distributions for the different signal modes with the background. In Appendix, we show distributions of some input variables for the signal and the background for the interested readers.

It should be remembered that relative importance or in other words the ranking of a variable might change for different production modes and also for different parameters like , , and so forth which can change the shape of the kinematic distributions. It is important to mention that this set of twelve variables used here may not be the optimal one. One can always improve the analysis with cleverer choices of variables.

In Figure 7, we show the BDT response by comparing two different production modes at a time. The and the fusion modes are very similar in nature and, therefore, it is extremely difficult to distinguish them. We do not consider the fusion further as it is very much identical to the mode. We show, by picking two production modes at a time, ten such possible BDT responses in Figure 7. These responses are substantially different for most of the combinations and therefore can be distinguished very efficiently. We observe that it is hard to distinguish the and the production modes as their BDT responses are not very different from each other. One should also notice that there are two peaks in the BDT response of the mode. This is because two types of different topologies, that is, the associated production and the vector boson fusion, contribute to the final state. In case of the bimodal distributions like this, one can use two different BDTs that are trained for two different topologies to further improve the analysis. This type of advanced analysis is beyond the scope of this paper.

As a side remark, one should always be careful about overtraining while using the BDT algorithm (or any other algorithm which uses nonlinear cuts). This can happen without the proper choices of the algorithm specific tuning parameters. One can check whether a test sample is overtrained or not by using the Kolmogorov-Smirnov (KS) statistics. Generally, if KS probability lies within the range 0.1 to 0.9 this guarantees that the test sample is not overtrained. For this purpose, one uses two statistically independent samples, one for training and the other for testing.

#### 4. Conclusions

Among the various resonance search channels at the LHC, the diphoton channel is particularly important as this channel provides a comparatively cleaner background. Generally, the diphoton resonance searches at the LHC assume that the resonance is produced from the fusion. Apart from the fusion production, many BSM theories predict TeV-scale scalars that decay to diphotons can dominantly be produced by other means, namely, through the quark-quark () fusion or through the gauge boson fusions (, , , and ). In this paper we consider an effective field theory of a heavy scalar that decays to diphotons. In this model-independent approach, the scalar can be produced in all the possible types mentioned above. We derive the exclusion limits on the mass and the effective couplings of the scalar using the latest 13 TeV ATLAS diphoton resonance search data with . While deriving the limits, we consider, for simplicity, only one effective coupling other than (since we only focus on the diphoton final state) which is nonzero. We have properly taken care of the modified cut efficiencies while recasting the limits set by the ATLAS collaboration. We find that when the scalar is dominantly produced from the fusion, the latest LHC diphoton resonance search data sets limit on the new physics scale TeV for the coupling for TeV.

In future, if a scalar resonance is seen at the LHC in the diphoton channel, the immediate important issue one has to investigate is how the scalar is produced. Some preliminary analyses have already been done in the context of the 750 GeV resonance where it is shown that the jet multiplicity distributions can be very different for the different production modes. In this paper we revisit the issue and show that the average jet multiplicity and the distribution can act as good discriminators. For better discrimination, we use a sophisticated multivariate analysis by combining twelve simple kinematic variables to distinguish one production mechanism from the other. Our analysis shows that one can identify different production mechanisms very efficiently at the LHC.

#### Appendix

#### Distributions for Signal and Background

For the interested readers, we show here the distributions of some input variables for the signal and the background. The signals in Figures 8 and 9 are for the and the fusion production modes, respectively. These distributions are obtained by applying the selection cuts defined in Section 3.3. The BDT responses for these two production modes with background are presented in Figure 10. We observe that the signal and the background distributions are very different in nature and, therefore, one could use a MVA to isolate the signal from the background. After filtering out the signal events from the background, one can use our method to identify the underlying production mechanism. It is expected that the signal distributions deviate more and more from the background as we increase the resonance mass. Therefore, isolation of the signal from the background becomes easier for heavier resonances. One can, therefore, tune MVA for lower masses and use the same optimized analysis for higher masses, for simplicity.

**(a)**

**(b)**

Notice that there is a second bump around 500 GeV and in the range 300–500 GeV in the background and distributions, respectively. Similarly, there is a second bump in the background distribution around 1000 GeV (this is expected since is correlated with the transverse momenta of the photons). This unusual shape of these distributions also leads to the bimodal nature of the background BDT responses in Figure 10. The origin of these peculiar second bumps in the background distributions is due to the selection cuts and used to obtained these plots. We have confirmed that these bumps go away with the removal of the above-mentioned correlated cuts.

#### Conflicts of Interest

The author declares that they have no conflicts of interest.

#### Acknowledgments

The author thanks Valery A. Khoze, Lucian A. Harland-Lang, Rikard Enberg, and Gunnar Ingelman for helpful discussions. This work is supported by the Swedish Research Council under Contracts 621-2011-5107 and 2015-04814 and by the Carl Trygger Foundation under Contract CTS-14:206.