The Di-Photon Excess in a Perturbative SUSY Model

We show that a 750 GeV di-photon excess as reported by the ATLAS and CMS experiments can be reproduced by the Minimal Dirac Gaugino Supersymmetric Standard Model (MDGSSM) without the need of any ad-hoc addition of new states. The scalar resonance is identified with the spin-0 partner of the Dirac bino. We perform a thorough analysis of constraints coming from the mixing of the scalar with the Higgs boson, the stability of the vacuum and the requirement of perturbativity of the couplings up to very high energy scales. We exhibit examples of regions of the parameter space that respect all the constraints while reproducing the excess. We point out how trilinear couplings that are expected to arise in supersymmetry-breaking mediation scenarios, but were ignored in the previous literature on the subject, play an important role.


Introduction
In the first presentation of LHC Run 2 data, both experiments ATLAS and CMS presented an excess in the di-photon mass spectrum for comparable invariant masses. The CMS analysis observed its largest excess in the di-photon mass spectrum based on 2.6 fb −1 of pp collisions at √ s = 13 TeV for an invariant mass of 760 GeV with a local significance of 2.6 σ and a global significance of smaller than 1.2 σ [1]. Similarly, the ATLAS collaboration reported the largest deviation from the background hypothesis for an invariant mass of 750 GeV using 3.2 fb −1 of data, leading to a local significance of 3.6 σ and a global significance of 2.0 σ taking into account the look-elsewhere-effect in the mass range of m γγ ∈ [200 − 2000] GeV [2].
After updating and refining their analysis, CMS achieved an improved sensitivity by more than 20 % and added a new data set which was taken with B = 0 T reaching as well a comparable 3.3 fb −1 . The modest excess at 750 GeV for the combined 8 and 13 TeV data remained with 3.4 σ (local) and 1.6 σ (global) significance [3]. ATLAS updated their 8 TeV analysis and confirmed the modest excess at 750 GeV in the Run I data set with a significance of 1.9 σ. Thus, the recent updates strengthen the hint for a new physics signal.
For the Spin-0 hypothesis and under the assumption of Γ/M S = 0.014 × 10 −2 (with M S the scalar singlet mass) the combined dataset of CMS with 3.3 fb −1 (13 TeV) and 19.7 fb −1 (8 TeV) gives the production cross-section times branching ratio into two photons to be σ 13 TeV · B γγ ≈ 3.7 ± 2fb. (1.1) while one analysis of the ATLAS data gives [4].
An interpretation of this excess is that it is due to the production and subsequent decay of a scalar resonance of mass 750 GeV; while there have been many alternatives proposed (too many to mention here), we shall restrict to that case here as the most obvious and least tuned option in perturbative theories. The existence of such a particle with a mass close to the electroweak scale implies a new hierarchy problem that cannot obviously have an anthropic explanation, and this naturally strengthens the case for low-energy supersymmetry. However, the observed rate of diphoton production via the resonance is too large compared to what is expected from a heavy Higgs companion of the light Standard Model (SM)-like one, and in particular it is very difficult to justify in the Minimal Supersymmetric Standard Model (MSSM) (see e.g. [5] 1 ). In fact, the interpretation of the excess is challenging for most previously proposed supersymmetric extensions of the Standard Model, and of the perturbative models proposed since the announcement almost all invoke additional vector-like fermions and/or bosons. For an early review see [8]. In this work we shall show, on the other hand, that a previously proposed supersymmetric extension of the Standard Model called the Minimal Dirac Gaugino Supersymmetric Standard Model (MDGSSM) [9] contains all of the ingredients to explain the excess.
Since the proposal in [10] of extending the MSSM with extra states in the adjoint representation of the Standard Model to allow Dirac gaugino masses, this possibility has been subject to many studies due to their theoretical and phenomenological advantages: they allow simpler models of supersymmetry-breaking due to preserving an Rsymmetry; their masses are supersoft [11] and supersafe from collider searches [12][13][14]; they ameliorate the SUSY flavour problem [15][16][17]; and contain new couplings which aid the naturalness of the Higgs mass [9,[18][19][20][21]. Indeed, multiple realisations have been proposed that differ by the fate of R-symmetry, the presence or absence of additional states and interactions [11, (for a short introduction see for example [58]). Here we consider the case of the MDGSSM which was introduced with a minimal content of extra states to automatically preserve unification of gauge couplings while allowing the new couplings to the Higgs to enhance naturalness and allow the boundary conditions to be unified at a high energy scale.
We will show that it is one of the most promising models when it comes to reproduce the diphoton excess. Without any ad-hoc addition, all the necessary ingredients are already present in the MDGSSM: • There is a singlet supermultiplet S introduced to give the Bino a Dirac gaugino mass. It is straightforward to identify its scalar (or pseudoscalar) component with the 750 GeV resonance.
• There are extra vector-like charged states, subsequently called "fake leptons" [59] as they carry the same quantum numbers as the Standard Model leptons. They were introduced in order to restore the automatic gauge coupling unification that was spoiled by the addition of the adjoint representations of the Standard Model gauge group. In this work, these states will increase the coupling of the scalar resonance to photons at one loop.
• There is an octet supermultiplet O required to give the gluino a Dirac mass. This contains colour-octet scalars which will generate a coupling of the singlet resonance to gluons at one loop (via trilinear scalar couplings), required for its production in gluon fusion.
One of the important constraints to impose on any new scalar S candidate to explain the excess is a bound on its mixing with the Standard Model Higgs. This mixing is not only induced at one-loop, but can be present already at tree level. The supersymmetric operator describing the Dirac gaugino bino mass leads to a modification of the U (1) Y D-term as where S R is the real part of S and ϕ j a scalar field with charge Y j under U (1) Y . Upon elimination of the auxilliary fields, this implies an interaction of the form: thereof a tree-level induced mixing. However, this is typically compensated by the presence in the superpotential of a term of the form: A precise evaluation of this mixing at the tree and one-loop level needs to be carried out carefully if one tries to identify the scalar S in models of Dirac gauginos with a 750 GeV resonance. Our parameter space is constrained by the requirements of stability of the vacuum avoiding existence of directions in the phase space of the model taking the fields expectation values to charge-and colour-breaking vacua. This is important as we shall see that trilinear terms will play an important role in generating the required amount of scalar production and decay into di-photons. Among the trilinear terms considered here some have not been explicitly discussed in the existing literature while they are expected to be generically present in the model. This is the case for example of soft terms mixing three adjoint scalars that we will show that they are generated in models of gauge mediation.
We shall keep couplings small enough to preserve perturbativity up to the GUT scale. This restriction can be of course relaxed if one allows for Landau poles below the GUT scale. However, as one of the virtues of the MSSM was to predict perturbative unification of gauge couplings, and was one of the motivations for introducing the MDGSSM, we shall place emphasis on finding regions of the parameter space which respect this condition.
To find the parameter space relevant for the diphoton excess we shall use the most sophisticated tool available: the code SARAH [60][61][62][63][64][65] and its SPheno [66,67] output. This is able to calculate the masses of all particles to full one-loop order, and twoloops in the gaugeless limit for the neutral (pseudo)scalars [68][69][70]. It can calculate renormalisation group equations of all couplings to two-loop order, including the masses and tadpoles in Dirac gaugino models as given in [71]. A guide to its use for studying the diphoton excess was described in [8]; we make some small modifications described in section 5.1. In particular, this will allow us to obtain the production and decays of our resonance at 8 and 13 TeV while simultaneously accurately computing its mass and assuring that the light Higgs mass is correct, and verifying that the mixing between the singlet and the Higgs is small (also computed at two loops). We shall find that quantum corrections to the spectrum of particles are not just important but essential for understanding how the model describes the excess.
Finally, we note that there have been three previous attempts to relate models with Dirac gaugino masses with the diphoton excess. In [72] as in this work the scalar component of S was the putative resonance; however, the entire coupling was driven by (1.4) which required very large Dirac gaugino masses (which would potentially flatten the Higgs potential). In [73] the candidate is a neutral component of a scalar doublet R 0 u introduced in the MRSSM to preserve R-symmetry, but the model required the R-symmetry to be broken to fit the excess and the Dirac nature of the gauginos played little role. As we were preparing to submit this work, [74] appeared, where the pseudoscalar component of S plays the role of the resonance; it couples entirely via superpotential couplings to coloured and charged fermions and thus requires large Majorana gaugino masses and charginos close to the threshold of 375 GeV to generate the couplings to photons and gluons. Here we will not require any Majorana masses, and will include only ingredients already allowed in the MDGSSM.
The paper is organised as follows. In section 2, we summarise the MDGSSM field content and interactions. To generate a large gluon coupling we require trilinear scalar adjoint couplings, the generation of which we describe in section 3 along with some observations on adjoint scalar masses. We discuss the constraints on the model in section 4; in particular, this includes a detailed study of vacuum stability, and an analysis of the constraints on colour octet scalars which are important and interesting in the context of this model. Our numerical results are provided in section 5 with some benchmark points to illustrate how our model reproduces the signal. Our results are summarised in the conclusions.

Model Content and Lagrangian
In this section we review the main ingredients of the Minimal Dirac Gaugino Supersymmetric Standard Model (MDGSSM) introduced in [9].

Field content
The MDGSSM field content can be seen as the minimal set providing the MSSM gauginos a Dirac masss while preserving two-loop unification and perturbativity of gauge couplings. We summarised it in Table 1. In addition to the chiral multiplets transforming under the adjoint representations of the gauge groups, it includes new fields charged under the lepton number global symmetry. They consist of extra Higgs-like doublets 2 R u , R d as well as two pairs of vector-like right-handed electron superfields E 1,2 in (1, 1) 1 andẼ 1,2 in (1, 1) −1 . Such states are compatible with an (SU (3)) 3 Grand Unification gauge group. This is the minimal set which enables a "natural" unification (unification without mass thresholds tuning) similar to the MSSM.
The adjoint chiral multiplets contain new complex adjoint scalars, S, T and O: I , T I , T −I , T +I are real scalars and pseudo-scalars, respectively.

Lagrangian
The superpotential for these fields can be written as where W Y ukawa contains the usual MSSM Yukawas part W DG contains the a priori R-symmetric contributions of the non-MSSM fields 3

Names
Spin 0 Higgs-like Leptons In this work we shall, as in [9], consider scenarios where R-symmetry is preserved by the superpotential (and thus these terms vanish). However we shall also consider the possibility that they do not vanish -so the superpotential violates R, in particular λ SO will play an important role in the following. For simplicity and to avoid lepton-flavour-violation constraints, we shall only the terms of the first three lines of (2.4) to appear with sizable couplings; the contributions of the last two must be small enough to be negligible for the purpose of this work, so we shall set them to zero throughout.
For the soft SUSY-breaking terms, from the MSSM we retain only the bilinear terms -i.e. conventional mass-squared terms and the B µ term. All the scalar triilinear and Majorana gaugino mass terms violate R-symmetry; while for B µ we suppose that, since R-symmetry is a chiral symmetry, we are breaking R-symmetry in the Higgs sector -and in fact it is only in combination with the superpotential terms mu, λ S , λ T R is violated. Hence in principle we can have an entirely R-preserving supersymmetrybreaking sector.
The soft SUSY breaking terms beyond those of the MSSM consist of 4 : • Dirac gaugino masses: (2.6) • soft terms associated with the adjoint scalars The terms on the last line have generally been neglected, but will play an important role in this work.
• soft terms involving the new vector-like leptons: Let us highlight that in an R-symmetry conserving model, one cannot simultaneously have the trilinears T SE (respectively T SR ) from (2.8) and the superpotential couplings λ SE (respectively λ SR ) from (2.4) as each term requires a different R-charge for the fieldsÊ andÊ (respectively R u and R d ) to be R-invariant.

Scalar mass matrix
We use the notationm where the effective masses for the real parts of S and T read: Then, at tree level the scalar mass matrix in the basis {h, H, S R , T 0 R } is [20]: where we have defined: which vanishes when λ S and λ T take their N = 2 values, (2.14) This matrix is diagonalised by the mixing matrix S ij . Of particular interest will be S 11 which measures if the lightest scalar eigenstate is Standard Model Higgs like, and S 13 which measures the proportion of the scalar singlet S R in this lightest eigenstate.

Generating Trilinear and Quartic Couplings
Previous studies of Dirac gaugino models have generally neglected the phenomenology of adjoint self-coupling terms, with an exception being a superpotential term κ 3 S 3 used in [19] to generate µ/B µ as in the NMSSM, and a recent brief discussion in [56]. In the case of superpotential terms such as λ SO these can be neglected when considering an R-symmetric visible sector; however, trilinear soft couplings such as T SO , T O (see (2.7)) are always allowed. It is therefore interesting to consider what values we expect from models of supersymmetry-breaking mediation.
Starting with a spurion analysis where supersymmetry is broken by either a Dterm D or F-term F , then if the mediating dynamics is at a scale M the terms in our effective Lagrangian should be given by powers of D M , F M , D M 2 , F M 2 with appropriate factors of couplings and κ l ≡ 1/16π 2 . Furthermore, quartic and higher-order couplings -which are "hard" SUSY-breaking parameters -are always generated, but do not lead to quadratic divergences because they appear suppressed by powers of the scale M which is the cutoff of our effective theory. Important in this work are the quartics such as L ⊃ λ 4S 24 S 4 which must have size λ 4S ∼ κ p l D M 2 q for some integer p, q (or similarly for F-terms with even q)); taking p = 1, q = 1 for a D-term we naively have a quadratic divergence in the scalar mass proportional to λ 4S but this yields In fact, this tells us that the case q = 1 is special because it implies a much larger correction at one loop than the direct mass, and could therefore destabilise the calculation. We shall return to this below.
As a first observation, if the mediation is by gravity, then M should be identified with the Planck scale (unless there is significant sequestering) and we should only consider the leading order terms. We would therefore require the quantum gravity theory to give us the terms T SO , T O at leading order D/M, F/M and the quartics must, by the above reasoning, be negligible.
On the other hand, in the case of low-scale supersymmetry breaking -where it was argued in [75] that this requires Dirac gauginos -√ F ∼ √ D ∼ M ∼ TeV, and we generate all terms at a similar order, which would include T SO , T O . However, the phenomenology is significantly changed by the presence of higher-dimensional operators and the goldstino couplings [50] and, since it is difficult to reconcile with perturbative unification, we shall not discuss this further here.
Finally, for gauge mediation M could be as small as √ F or √ D but there is no a priori upper limit on M until we choose a particular quantum gravity embedding. The Dirac gaugino masses are expected to be generated at one loop and be of order κ l D M or κ l

Adjoint couplings in gauge mediation
One of the most interesting issues in the construction of gauge mediation models with Dirac gaugino masses has been that of the adjoint scalar masses: in the simplest realisation, only a B-type mass-squared L ⊃ − 1 2 B Σ Σ 2 is generated at leading order in D/M 2 , and not a conventional mass-squared L ⊃ −m 2 Σ |Σ| 2 . This happens for one pair of vector-like messengers Q,Q having charges under a hidden U (1) of +1, −1, where the U (1) obtains a D-term. This was noticed from the earliest models [25,76] with the original proposed solution being to add a supersymmetric mass for the adjoint -which would also violate the R-symmetry and generate Majorana masses for the gauginos, with a see-saw effect. However, an alternative solution was found to be to introduce additional messenger states with non-diagonal couplings to either the adjoints (in the D-term case) [28,31] or an F-term spurion [26,28,31]; in the D-term case this requires the couplings to violate the U (1)-charges. In [31] examples were given where the ratio of B-type to conventional masses is arbitrary. The general ansatz was to couple the adjoint to messenger fields Q i ,Q j and to possible F-term spurions X via superpotential couplings and D-terms via charges e i ,ẽ i which we can write as a matrix e ij (Q i Q * j −Q iQ * j ). More recently, the issue has been re-examined. One suggested approach, dubbed "Goldstone gauginos," is to promote the adjoints to be the Goldstone bosons of a broken symmetry [54,55]; however, this solution would lead effectively to no higherorder interactions for our adjoint scalars and we do not consider it here. More in the spirit of the earlier works, the issue was rephrased in the language of effective operators in [33,47,56], where it was claimed that the explanation for the absence of conventional mass-squared terms for the adjoints at leading order is that the operator responsible for the generation of a leading-order mass-squared term should be where ψ,ψ are a pair of fields charged under the hidden U (1) with charges ±q which obtain vevs (and thus generate a contribution to the hidden D-term). The above operator is generated by including terms in the superpotential that mix the messengers Q,Q with other pairs of fields N,Ñ which are neutral (or at least have different charges) under the hidden U (1), so that the vevs of ψ,ψ generate messenger mixing terms. This is clearly nearly equivalent to the above ansatz, and can be written in the form where we now write the mass terms as violating the U (1) charges instead. If we start with the case of no couplings/mass mixing terms that violate the U (1) D-term charges, we shall first give a simple proof that the conventional mass term |Σ| 2 vanishes at leading order for any number of messengers, and then look at higher-order terms. Considering first the visible gauge group to be U (1), we have the effective potential contribution from the messenger scalars (since the fermion potential is independent of D): are the supersymmetric mass-squared matrices. Now, if we take the couplings to preserve the U (1) charges then we can write the since the eigenvalues of M 2 Q and M 2Q are equal. Next, by taking the derivative with respect to Σ we find only a holomorphic function of Σ: As an example, consider the simple model of a single messenger where the matrices become numbers; then we havẽ This potential manifestly has trilinear and quartic couplings, although at order respectively. Indeed, if we continue with the ansatz (3.1) then it is easy to see that there are no terms of linear order in D, because M 2 Hence to have large cubic interactions we should start from ansatz (3.3). In this way, in order to have an interesting phenomenology we require either D ∼ M 2 with both at a low scale, or we require (as proposed in [47]) that with some cancellation between the two terms so that we can have m Σ ∼ T SO . Note that once we take this ansatz with non-zero mixing between the messengers and [M, e] = 0 we typically generate trilinear terms in the potential -but also tadpoles. The issue of tadpoles is then easily circumvented by embedding the coupling of the singlet adjoint S to the SU (3) and SU (2) adjoints into the generator . This then also means that the couplings of the singlet adjoint S are related to those of T and O; for example, for T SO , if we have calculated the coupling for U (1) messengers as being L ⊃ 1 6 T Σ Σ 3 , then we have where T a 3 , T b 4 Constraining the MDGSSM from the diphoton excess We analyse in this section various theoretical and experimental constraints lying on the general model presented above. We start by considering the basics of production and decay of the scalar singlet and then study the most relevant collider constraints on our model. We finally investigate the requirements we need to impose in order to remain perturbative up to the GUT scale and avoid the appearance of Charge or Colour Breaking Vacuua.

Production and Decay in the MDGSSM
In the narrow width approximation in which the mediating S singlet is automatically on-shell, we can approximate the cross section of the complete process pp → S → γγ as follows: Assuming a spin-zero particle produced resonantly via gluon fusion, we arrive at taking C 8TeV gg = 174 and C 13TeV gg = 2137 as values arising from the parton distribution functions [77], respectively. An important aspect of our calculation is that for a more realistic estimation, we have taken into account the K-factors K 8,13 for the full N n LO production of H + jet compared to the tree-level process. We have estimated K 8 1.9 from the comparison of the leading-order effective vertex from MadGraph and the Higgs Cross-section working group value for a Standard-Model-like Higgs of 750 GeV at 8 TeV. We will take conservatively the same value for K 13 .
Let us first consider the coupling to two gluons. The process S → gg is a priori generated by loops of squarks, scalar octet and gluinos. The amplitude is of the form , the sums runs over all scalars and fermions, and .
and Q f , Q φ , g Sf f and g Sφφ are the electric charge and coupling with the singlet of the fermions and scalars participating in the triangular loops. The loop fonctions A S 0 and A S 1/2 have a maximum at the resonant mass M S /2 ∼ 375 GeV. We will therefore generically require masses close to this scale in order to enhance the cross-section. The main contributions to the loop will be: • D-term-induced couplings between the squarks and the singlet, generated by the Dirac masses operator of Eq. (2.6). Theses couplings are proportional to the hypercharge of the squarks and the Dirac mass m 1D . They are sizeable only for large Dirac mass m 1D .
• Soft terms trilinears couplings from (2.8) between the adjoint scalar octet and the singlet. They give a sizeable contribution but unfortunately are strongly constrained from vacuum stability bounds.
A priori, one could have expected a contribution from the Dirac gluinos. However, we observed that pure Dirac gluinos do not contribute at all to the amplitude. This remark is of crucial importance for the pseudo-scalar S I which can only couple to gluons through fermions loops as we assume CP-conserving interactions. If no Majorana masses for the original gluinos are introduced, the pseudo-scalar is practically not produced. However, if we allow for the presence of an additional Majorana mass term, the pseudo-scalar S I can then participates in the S → γγ cross-section, potentially leading to a "double-peaks" scenario, as we will see later.
We now turn to the amplitude to diphoton. This is given for a scalar by In order to get an idea of the enhancement we need from the square term, let us find the smallest value of Γ γγ leading to a σ(S → γγ) 2 fb. In the limit in which Γ gg dominates the decay width, we can use Eq. (4.2) to get which is an order of magnitude bigger than the numerical factor in (4.4). The key issue will therefore be to populate the sums in the square terms of (4.4) since the amplitude will very roughly scale as N 2 , with N the number of particles participating in the loop. The main contributions will come from • D-term-induced couplings between the sleptons and the singlet, they are again proportional to the hypercharge of the sleptons and to the Dirac mass m 1D . They are therefore sizeable only for large Dirac mass.
• Superpotential-induced couplings between the fake leptons and the singlet from the terms of (2.4) in section 2. They are the main contributions in our model. 5 • Soft terms trilinears couplings from (2.8) between the fake sleptons and the singlet. They are again strongly constrained from vacuum stability bounds.
An important remark here is that the two last contributions are mutually incompatible in presence of a preserved R-symmetry as we already stressed in Section 2.1.

Constraints from Higgs mass mixing and 8 TeV data
A crucial property of the singlet S is that it will in general mix with the Higgs eigenstates. This is in our case an undesirable feature since it will lead to tree-level decays of S into tops, W , Z or Higgs which could easily overcome the one-loop decay into photons.

Analytical Estimate
Building on the notations introduced in the previous sections, we can use the minimisation condition of v S on the off-diagonal element ∆ hS of the scalar mass matrix given in (2.13) to find (see [78]) where we used the effective mass parameter From this basic analytical calculation, we see that we can minimise the tree-level mixing by choosing: In general, this relation will be modified at one-loop, but the property that one value of λ S is favored will remain and is easily observable in our coming Figures.

Experimental Bounds and Naturalness
Such a mixing with the Standard Model Higgs will modify the Higgs sector observables. From [79] we find the latest constraint on the 125 GeV Higgs global signal strength µ average to be µ average = 1.09 +0.11 −0.10 , (4.9) In our case this is modified by a factor of |S 11 | 2 , where S is the mixing matrices of the scalar sector; the above constraint gives us This condition is in fact satisfied quite easily, as can be seen from Figure 1 where we show the contours for the Higgs mass and the mixing matrix element S 31 as a function of tan β and λ S . An important comment regarding this Figure is that a 125 GeV also favors small mixing.  The thin black lines represent the 2% and 4% mixing contour lines. The anomalies around tan β ∼ 2.5 corresponds to the region where the two-loop effective potential used to determined the Higgs mass suffers from the so-called "Goldstone boson catastrophy" (see [69] for more details).
More stringent constraints arises from the non-observation of any excess in the 8 TeV data for the ZZ, and hh, dijets and W W channels. As the mixing between S and h induces a tree-level decay one naively expect a percent-level suppression to be necessary. Since the S is mostly produced by gluons fusions in our scenarios, we request that (see [77]): which gives the most stringent constraints on the mixing between S and h. The di-Higgs channel is proportional to the tree-level mixing term without passing through the mixing, because the vertex is given by ∆ hs /v (plus smaller terms proportional to the mixing matrix elements S 13 , S 31 ); we have and Neglecting v T and mixing with the triplet as small effects, we can then write Translating these into constraints, we see that it is the Z decays which are most important.
Notice that the only loop decays included in this paper are S → γγ and S → gg (as they do not have a tree-level contribution). A priori in the negligible mixing region, one should also consider the other diboson loop decays (in particular to Zγ). However, almost all of the new fields contributing to the loop decays will be SU (2) singlets so that the decay to diphoton will be the dominant diboson decay channel. The only exceptions are the new doublets R u and R d which should mostly decay to W W , ZZ and Zγ. Due to the interference with the tree-level processes the loop contribution to these processes is not currently calculated in SARAH; their implementation is eagerly awaited in future work, but here we note that they will not have a significant impact on our results as described in [8].
Finally, the VEV of T gives a contribution to the W boson mass and the electroweak precision data give bounds on it. One must examine the induced correction ∆ρ to the Veltman ρ-parameter: with ∆ρ given analytically at tree-level by ( [78]) where v is the usual Standard Model Higgs VEV. In order to be below the experimental constraints, we need ∆ρ (4.2 ± 2.7) × 10 −4 , ( [78] -see also [21,49] withm 2 T R = m 2 T + 4m 2 2D + B T , therefore, small ∆ρ require large triplet Dirac and soft masses. This requirement can often be at odd with naturalness which prefers smaller triplet masses. Indeed, radiative corrections induced by the adjoint triplet scalars to m 2 H u,d are [78]: with Λ the UV cut-off, m 2 H u,d , m 2 T the squared masses for Higgses and scalar triplet T , and λ T the coupling defined in (2.4). For Λ at the Planck scale, requiring a fine-tuning ∆ T = δm 2 H /m 2 H better than 10% finally gives us In Figure 2, we show the allowed region for λ T and m 2D for m T = 450 GeV. ∆ρ has been obtained at one-loop using the Spheno [80,81] code generated by SARAH (see ref. [82][83][84][85][86]). We see that the Higgs mass prefer large values of λ T but that the following three requirements are perfectly compatible: (1) a 125 GeV Higgs, (2) a natural mass for the triplet and (3) a parameter ∆ρ smaller than the current constraints.

Bounds on colour octets
In this work we shall be interested in the case when either the scalar or pseudoscalar colour octets are lighter than a TeV. Even though such light scalars should be copiously  produced in pairs at both 8 and 13 TeV, as shown in figure 3, their decays are loop suppressed and this inhibits single production.
Since current limits place all squarks above about 800 GeV, then, as first discussed in [27,87], the octets decay only to gluons and quarks -in particular almost entirely top quarks. This means that the possible signatures are four jets, dijet/ditop searches, and four tops. Up until relatively recently the constraints on them were rather weak, with dijets providing no constraint, and a mild constraint from ditops [88]. However, now the four top channel is particularly important: [89] placed a limit of 32 fb at 8 TeV, and [90] found 370 fb for Standard-Model-like kinematics, or 140 fb with and EFT pointlike interaction, at 13 TeV.
To interpret the implications of these searches for our model, we could in principle do a full recasting along the lines of [91]; however, for simplicity we shall consider instead the cross-section times branching ratio approach, taking the most conservative values of twice the tree-level cross-section (i.e. a K-factor of 2) and a limit at 13 TeV of 140 fb. To compute the branching ratio into four tops, we require the widths into gluons and tops; while expressions were given for these originally in [27,87], those papers used complex octets, which is not appropriate for our case where the necessarily large ( 2 TeV) gluino mass causes a large splitting. Instead we require the expressions presented in [50], which we shall not reproduce here but to which we refer the reader.
The first important observation is that the pseudoscalar octet does not couple to gluons, and so pair production of pseudoscalar octets yields only four-top events, and by our above criteria excludes pseudoscalars below about 880 GeV by the 13 TeV data. These are therefore less interesting for our analysis.
On the other hand, the scalar octet couples to squarks via its D-term coupling, and so couples to gluons. Since it couples to all coloured squarks, this can potentially be large. However, to be very conservative, we show production times branching ratio of four-tops via scalar octets in figure 4 at 8 and 13 TeV with the limits shown using a K-factor of 2, as we vary the octet mass and for three different values of the Dirac gluino mass, where the first two generations of squarks are decoupled (i.e. heavy and degenerate). To produce these, we take left-handed stops and sbottoms of 1200 GeV, right-handed stops of 800 GeV, and decoupled right-handed sbottoms (at 4 TeV). We neglect all squark mixing (which is a good approximation in this model). Since the couplings involve a cancellation between left-and right-handed squarks, this is very conservative: if we took heavier left-handed squarks, we would enhance the gluon rate relative to the top rate (because it has a contribution from sbottoms as well as stops) weakening the bounds.
We conclude that for 2.5 GeV gluinos, the octet scalars must be heavier than 500 GeV; but for 3 TeV gluinos there is no constraint.

Perturbativity and Landau Poles
The field content of the MDGSSM, and more precisely the two pairs of vector-like electronsÊ andÊ as well as the doublet R u , R d , have been chosen to have one-loop unification by completing the 8 0 + 3 0 + 1 0 set of adjoint multiplets into a complete GUT representation of (SU (3)) 3 (see [9]). We have furthermore checked numerically that gauge couplings remain safely perturbative at two-loops up to the GUT scale, consistently with the results of [9]. Once the GUT scale is determined, we require perturbation theory to be valid up to the GUT scale. We choose as perturbativity requirement that all Yukawa couplings should remain smaller than √ 4π. As we will see now, this gives strong constraints on the Yukawa couplings. At one-loop, the beta functions for λ SE , λ SR , λ SO , λ S and λ T form a coupled system given by: where the dots contain the contributions from the other couplings. Before studying this system numerically, we point out some peculiarities of these expressions: • The gauge couplings contribute negatively to the beta function, increasing the stability. In particular, λ SO is strongly stabilised.
• In the limit λ S → 0, λ T completely decouples from the other Yukawa couplings.
• The perturbativity of the coupling λ S will be critical as: (1) the gauge couplings and top Yukawa already give a positive contribution ∼ 1.1 to its beta function; (2) all the other Yukawas feed intro its beta function and conversely λ S feeds into all the beta functions.
We have numerically constrained the initial values for λ SE , λ SR , λ SO , λ S and λ T at the low scale (SUSY scale), so that they remain perturbative up to the GUT scale. We use the two-loop RGEs generated by the public code SARAH (see ref. [82][83][84][85][86] and ref. [78]).
In Figure 5, we study the case of λ SO = 0, which will be relevant for the two R-conserving scenarios R a and R b . The perturbativity bounds are shown in the planes λ S /λ SE and λ S /λ T . As expected, we obtain the strongest constraints for λ S , especially in the large λ SE case, which is the one of interest in this paper. Furthermore, we recover that for λ S → 0, λ T is insensitive to the other Yukawa couplings. Adding the parameter λ SO further constrains the Yukawa couplings. This is shown in Figure 6 where we present the perturbativity bounds on λ SE and λ SO for various values of λ S and λ T . We see that for λ SO ∼ 0.65, one should take λ SE < 0.65 to be safely perturbative. Furthermore, as expected from the one-loop beta functions, λ SO has an increased stability thanks to the strong gauge coupling contribution, allowing values up to 1.4 for low λ SE . Notice in the right-hand plot of Figure 6 that in the limit λ S → 0, we recover that λ T decouples from the other Yukawas.

Vacuum stability
We now turn to the constraints from vacuum stability; since we have significant trilinear scalar couplings then this is of crucial importance. The tree-level scalar potential can be decomposed into four main contributions: = . = .
= . Figure 5: Perturbativity bounds on our model, around the first benchmark point from Table 2, obtained from the requirement that no couplings overtake √ 4π before the GUT scale. We consider λ SR = λ SE . Left plot: Bounds for (from left to right) λ SE = 0.7, 0.5, 0.3, 0.1 in the λ S /λ T plane, all points above the curves are excluded. Right plot: Bounds for (from left to right) λ T = 0.9, 0.7, 0.4, 0.1 in the λ S /λ SE plane, all points above the curves are excluded.   Table 2, obtained from the requirement that no couplings overtake √ 4π before the GUT scale. We consider λ SR = λ SE . Left plot: in the λ SO /λ SE plane with from left to right (λ S = 0.3, λ T = 0.7) and (λ S = 0.05, λ T = 0.85); all points above the curves are excluded. Right plot: in the λ SO /λ T plane with from left to right (λ S = 0.3, λ SE = 0.65) and (λ S = 0.05, λ SE = 0.65); all points above the curves are excluded.
with V g , containing the D-term contributions, V W the superpotential contributions and V soft the soft SUSY-breaking terms. The final term V hard consists of "hard" dimensionless quartic terms that are generated at the SUSY-breaking scale and look like hard SUSY-breaking terms discussed in section 3.
We have where ϕ j are the scalar components of the matter chiral superfields, possibly in the adjoint representation and M a j is the matrix of the gauge representation of ϕ j . Let us leave aside the triplet contribution (we are considering a heavy triplet and therefore expect a near-zero VEV for it) and focus on the singlet and octet terms. Similarly, we will leave aside the squarks contribution as we are not considering large A terms and therefore do not expect them to acquire a color-breaking VEV. We have then with (T a ) bc = (−if abc ) and f abc the SU (3) structure constants. We now turn to the superpotential contributions (we suppress the i indices forÊ i andÊ j and the "·" denotes SU (2) indices contraction by tensors) and find: The only "hard" SUSY-breaking terms that will be of relevance to us will be a quartic octet coupling: which is of course not the only such possible term but is the most important. After adding the soft and hard SUSY-breaking terms, we obtain with 24) and the mixed contributions

Charge-breaking minima
First, we investigate the S,Ê,Ê sector, which can drive charge-breaking minima when they all acquire a vev. The relevant tadpoles (for just one pair ofÊ,Ê ) are We have two limits relevant for this model: relevant for this paper is the case that m DY is not large, in which case the most dangerous direction is the "classic" D-flat direction |Ê| 2 = |Ê | 2 . When S,Ê,Ê develop vevs, we can decompose the complex fields into real and imaginary parts; without loss of generality we can putÊ =Ê ≡ Solving then the equation for the singlet tadpole, we find the where we definedT Clearly we observe that we have appearance of a charge-breaking vacuum if However, it is only lower than our vacuum if the weaker condition is satisfied, or equivalently The analagous constraints also apply for the pseudoscalar direction, and also for the S, R u , R d sector.

Colour-breaking minima
A crucial part of our analysis is the presence of trilinear couplings of the singlet to the octet, which generate a coupling to gluons. However, just as the couplings to the selectron-like states allow charge-breaking minima, the octet scalar couplings permit colour-breaking minima. The analysis is identical for O R or O I with the opposite sign for T SO , so let us choose O R . The tadpole equations read We therefore see that the supersymmetric terms are equivalent to putting λ O = λ 2 SO , λ H SO = λ 2 SO ; in an analysis identical to the previous subsection we find that an additional vacuum exists when but that the minimum is only lower than the colour-preserving one when 33) or equivalently, when λ SO = 0, If we choose to break R-symmetry only in the Higgs sector via a B µ -term, then λ SO = 0. In this case, we need to rely on λ H SO and λ O only to stabilise the potential, leading to very strong constraints on the trilinear T SO . For instance, if we have λ H SO , λ O ∼ O(0.04) and a 400 GeV scalar octet, the trilinear T SO must be smaller than 310 GeV to ensure that the colour-preserving vacuum is stable.

Prelude
While the MDGSSM has a large set of free parameters, the most relevant ones can be divided into three roughly independent sets controlling different features: 1. Higgs and singlet masses and mixing: m 1D , m S , B S , tan β, µ, λ S and λ T . 3. Singlet decay amplitude to γγ: T SE , T SR , λ SR , λ SE supplemented with soft masses and B terms for the fieldsÊ,Ê , R u and R d .
The first set is dedicated to reproducing the measured Higgs boson mass as well as a 750 GeV scalar singlet. The value of λ S need to be adjusted to have a small mixing between both scalars, which is necessary both for the diphoton cross-section and for having m H ∈ [122, 128]. The second set can is then used to enhance the production rate of singlet through gluon fusion. The trilinears T SO are crucial in this respect as they allow the scalar octet to participate in the loop-induced coupling Sgg, greatly increasing the singlet production rate. Finally, the last set of parameters is used to increase the diphoton amplitude. The superpotential Yukawa couplings λ SE and λ SR from (2.4) are constrained to be below 0.7 to avoid the appearance of Landau poles before the GUT scale. The trilinears are mainly constrained by enforcing that the scalar fieldsÊ,Ê , R u and R d does not get a charge-breaking vacuum expectation value. We will investigate various scenarios that we can classify according to the presence or not of the R-violating terms (2.5): -R a : where we will consider large Dirac mass m 1D , so that the coupling to gluons proceeds through squarks loops. -R b : where instead consider small Dirac masses but light scalar octet, so that the coupling to gluons proceeds through scalar octet loops.
• R-symmetry violating models, for which we can have additionally the terms (2.5) and the trilinears T SE and T SR . We consider -/ R a : A generalisation of scenario R a with λ SO and the trilinears T SE and T SR included. -/ R b : Similar to the model / R a , but we further tolerate the presence of a Majorana gauginos mass terms. This allows to simultaneously produce the scalar S R and pseudo-scalar S I singlet and have a "double-peaks" resonance set-up.
In the following we shall present results of a numerical investigation of the parameter space of the MDGSSM for various scenarios. To do this we used the package SARAH to produce SPheno code to calculate the spectrum, production rate and decays. We created a new model file for the MDGSSM including the adjoint couplings λ SO , T SO . However, we found that modifications to the SPheno code were necessary: 1. We use pole masses instead of DR-masses for the selectrons and octet scalars in the calculation of loop couplings with the neutral scalars. This is because these masses are the most important for the gluon and photon couplings of our 750 GeV candidate, and can differ by more than a factor of two; as described in [8], using the DR masses is less accurate and so we employ pole masses just for these particles.
2. To facilitate our search for valid parameter points, we produced two different versions of the code. The first solves the tadpole equations for mass-squared parameters m 2 Hu , m 2 H d , m 2 S , m 2 T taking v S , v T as inputs; while this is the appropriate choice for implementing the loop corrections to the scalar masses correctly, it is, however, difficult to choose the vacuum expectation values v S , v T (since loop corrections can rapidly change the values of m 2 S , m 2 T by several orders of magnitude). We therefore use this version of the code to check the results of our second code, which was specially modified to first solve the two-loop tadpole equations numerically for v S , v T , and then compute the tadpoles and masses using these values as inputs, solving for m 2 S , m 2 T again along the same lines as the first code. While this is computationally expensive (computing the two-loop corrections to the neutral scalar tadpoles twice for each point) it is the most efficient way to correctly identify points -and not miss points where, for example, m 2 S may be identified as tachyonic at "tree level."

R-Symmetry conserving Scenarios
Consider first the scenarios R a and R b where we include only the R-symmetry conserving adjoint couplings. Under these constraints, the singlet production proceeds mainly by gluons fusion through loops of squarks (controlled by g Y m 1D ) and (pseudo-)scalar octets (controlled by the trilinear T SO ).

Squark-induced gluon fusion
We start with scenario R a and present in Table 2   The main aspects of this scenario are the following: • We limit the R-symmetry breaking to the Higgs sector, and therefore choose Rcharge of the fake fields to allow λ SE , λ SR superpotential terms but not trilinears T SE /T SR and the corresponding B-terms.
• The loop coupling to gluons will proceed through squark loops, with singlet/squarks coupling enhanced by a large Dirac mass m 1D • The loop coupling to photons have numerous contributions through loops of fake fields (both fermions and scalars) and sleptons.
• Finally, because of the large Dirac mass, one need a sizeable negative B S to ensure that the scalar singlet has a mass of 750 GeV.
Notice that a satisfying feature of this scenario is that we do not need to fine-tune the mass of the fields participating in the loop coupling beween S and γγ.
Regarding the scalar singlet production, gluon fusion proceeds mainly through loops of 800 GeV right-handed stop and TeV right-handed first two squarks generations, while left-handed squarks are heavier at 1.75 TeV. As a consequence, the mass of the right-handed stop is a critical parameter in enhancing σ γγ , we illustrate this dependence in Figure 7 where we plot the S → γγ cross-section as a function of the stop one-loop mass, by varying around the benchmark point of Table 2. We see that the cross-section decreases very rapidly with the stop mass.  Regarding the scalar singlet decay to diphoton, it proceeds both through loops of light right-handed sleptons (we consider left-handed sleptons above the TeV) controlled by g Y m 1D and loops of fake leptons,Ê,Ê , R u and R d which are controlled by a unified Yukawa λ SR = λ SE = 0.7. Furthermore, the fake sleptons also contribute with couplings controlled by g Y m 1D . In order to maximise the overall contribution, one has to take care that no cancellations occur between the various contributions (particularly for the D-term-induced couplings, which are proportional to the hypercharge of the scalar participating in the loop). Refering to Table 1 we see that one possible choice is lightÊ, R d and right-handed sleptons and heavierÊ , R u and left-handed sleptons.
In order to have sizable contributions from the (fake) sleptons, we need a reasonably large singlet Dirac mass m 1D ∼ 1250 GeV, this has the added benefit that it also enhances the squark contributions to the scalar singlet production rate. On the other hand, it increases the tuning of λ S necessary to have a small mixing and additionally implies that we have either a small tan β or a somewhat large µ term as can be seen from Eq. (4.6).
Overall, Figure 8 presents the cross-section obtained in the λ S /µ E plane by varying around the benchmark point of Table 2. Roughly speaking, this figure combines on the abscissa the constraints from mixing with on the ordinate the requirement that the particles participating in the loop have masses close to half that of the resonance. 6 We see from Figure 8 that the main requirement in our model is that we must consider values of λ S tuned at the level of a few percent. We can see that the constraint from the ratio Γ ZZ /Γ γγ is significantly weaker than the requirement on the crosssection.   Table 2. The black contour shows the most constraining ratio from (4.11) while the red contours shows the pole mass for the fake leptons.

Octet-induced gluon fusion
Let us now consider the case R b where we take a light scalar octet. Following the discussion of the previous section, its mass is not constrained as long as we take a large Dirac gluino mass. We will therefore focus on m 3D = 3 TeV. Since m 3D contributes at tree-level in the mass of the scalar octet, we require large negative B O ∼ −4m 3D in order to have it close to the resonant mass of 375 GeV. While this is a new source of tuning, the fact that the scalar octet provides a sufficient coupling between the gluons and singlet means that we no longer need a sizeable Dirac mass m 1D as in R a . As a consequence, the tuning on λ S is milder in this scenario, as can be seen from Figure 9. We have presented in Table 2 a benchmark point for this scenario.
As the singlet Dirac mass is small, the sleptons do not contribute to the singlet decay to diphotons, in stark contrast with scenario R a . One relies on loops of fake (s)leptons to increase σ S→γγ . The crucial parameter in this model is therefore the fake leptons mass, as we illustrate in Figure 9.   Table 2. The black contour shows the most constraining ratio from (4.11) while the red contours shows the pole mass for the fake leptons.

R-Symmetry violating Scenarios
If we do not constrain ourselves to an R-symmetric scenario, further regions in parameter space open up which can lead to an enhanced di-photon signal. As already stated in Sec. 5.1, we discuss mainly two different R-violating scenarios. In the first scenario, we consider the trilinear coupling κ still to be small, but allow for the trilinear couplings T SE and T SR as well as for the Yukawa-coupling λ SO . With the presence of the latter (λ SO = 0.65), it allows us to further increase the trilinear coupling T SO (e.g. from 200 GeV in Scenario R a to 1.5 TeV in Scenario / R a ) without leading to colour-breaking minima (cf. Sec. 4.5). As can be seen in Fig. 10 partial decay width of Γ S→gg , as it increases the coupling of the pseudo-scalar octets (m σ 0 = 886.30 GeV) in the loop to the singlet. However, this effect is reduced by an increase of the octet mass via loop effects. At the same time, the increase of T SO leads as well to an decrease of the mass of the fake sleptons, which increases the partial decay width of Γ S→γγ . The production via squarks (mq > 1.7TeV) is suppressed in this scenario. A further enhancement of Γ S→γγ is achieved by increasing T SE = T SR . The fake sleptons and sfermions are below 400GeV thus leading as well to an enhanced partial decay S → γγ. Similar to the R conserving Scenario, we account for a mass hierarchy betweenÊ, R d andẽ c andÊ , R u and L to prevent cancellations from D-term induced couplings, which could still be further increased to allow for an even larger production cross section. Table 3 indicates further input parameters and shows that this benchmark point is in full agreement with current experimental exclusion limits and features a production cross section of σ(S → γγ) = 3.1 fb. Choosing a large m 3D = 1600 GeV, the gluino mass lies above 1800 GeV, and m T = 1250 GeV guarantees that the ρ-parameter is below the exclusion limits. We have further checked that such a scenario is not yet excluded by resonant production of W W , ZZ, hh or gg. : λ S − µ E plane for benchmark scenario / R a . Fig. 11 shows the λ S − µ E plane around the benchmark scenario, where a similar behaviour as for the R-conserving scenarios can be observed. Although one would naively expect that in the R-violating scenario the necessary tuning would be less pronounced, the scenario is further constrained when requiring perturbativity of all Yukawa couplings up to the GUT scale. As discussed in Sec. 4.4, including λ SO as free parameter, it further constrains the maximal value of the other Yukawa-couplings. In the studied scenario, a Yukawa-coupling of λ SO ≈ 0.65 implies an upper bound on λ SE < 0.65 as well. However, generally, larger values of λ SO (and respectively of λ SE ) are possible.
In the second R violating scenario, we want to discuss the possibility of having a degenerate spectrum of the scalar and pseudo-scalar singlet. Only when the soft SUSY-breaking mass term M 3 is included, production via gluon fusion is possible and leads to a sizeable mass splitting between the Majorana gluinos. Having two particles leading to a di-photon signal, a broad parameter space opens up. This can for example be seen in Fig. 12 showing the corresponding λ S − µ E plane around the corresponding benchmark scenario / R b . With the pseudo-scalar singlet not being constrained by mixing with the Higgs like for the scalar singlet, a large variation in λ S is possible as with respect to the R conserving scenarios or / R a . Due to this degenerate scenario with the pseudo-scalar being less constrained, the requirement of having lower pseudo-scalar octet masses for boosting the production via gluons is loosened (m σ 0 = 886.3GeV), as well as for the masses for the fake sleptons and selectrons.

Conclusions
The MDGSSM is promising as a model that reproduces the di-photon excess observed at both LHC experiments, ATLAS and CMS. It automatically contains a singlet with both scalar (S R ) and pseudoscalar (S I ) components that can both be at the origin of the resonance. It is quite easy to fix the mass of one or both of them at 750 GeV. In the latter case, a small splitting will simulate the larger width of the resonance for which there is a mild preference in the present ATLAS data. Also the model contains new states beyond those of the MSSM, triplets, octets and fake leptons, that can be used in the loops to generate both the production of the singlet and its decay to photons. We have shown that there are diverse experimental constraints that are quite stringent. We have found that if the resonance is to be identified with the scalar S R , keeping its mixing with the Standard Model Higgs within the experimentally allowed range represents the most constraining issue. We have found that a certain amount of cancellation is needed between certain parameters and this can be translated in a tuning of the trilinear λ S at the level of a few percent. Fortunately, this happens to values of λ S sitting in a quite natural range, near the values expected from an N = 2 supersymmetric origin of the coupling.
We found that while remaining within the assumptions of the MDGSSM -perturbative couplings up to the GUT scale, R-symmetry-breaking only in the superpotential -the signal can be easily fit by including new dimensionful trilinear couplings of just the adjoints. The latter have not attracted attention in existing literature on Dirac gaugino models, despite the fact that they respect R-symmetry and so are always allowed. While they come out typically small in some scenarios of supersymmetry breaking, this is not always the case and they are expected to be present in the model. We have provided a first comprehensive discussion on this point. We then performed numerical scans of large parts of the parameter space using the most advanced tools available and in particular the most sophisticated calculation of the Higgs mass (up to two loop order). We found different regions of the parameter space of the MDGSSM (and given example benchmark points) satisfying all existing constraints while providing a good fit to the observed di-photon excess. Moreover, the relevant points have large quantum corrections (in particular to the singlet mass and vacuum expectation value) underlining the importance of using numerical tools.  Table 3: Overview of the input parameters, physical masses, constraints and production cross section for both R-violating scenarios / R a and / R a .