Hard Photons in Hadroproduction of Top Quarks with Realistic Final States

We present a complete description of top quark pair production in association with a hard photon in the dilepton channel. Our calculation is accurate to NLO in QCD. It is based on matrix elements for $e^+\nu_e \mu^-\bar{\nu}_\mu b \bar{b}\gamma$ production and includes all resonant and non-resonant diagrams, interferences, and off-shell effects of the top quarks and the $W$ gauge bosons. This calculation constitutes the first full computation for top quark pair production with a final state photon in hadronic collisions at NLO in QCD. Numerical results for total and differential cross sections are presented for the LHC at a centre-of-mass energy of $\sqrt{s}=13$ TeV. For a few observables relevant for new physics searches, beyond some kinematic bounds, we observe shape distortions of more than $100\%$. In addition, we confirm that the size of the top quark off-shell effects for the total cross section is consistent with the expected uncertainties of the narrow width approximation. Results presented here are not only relevant for beyond the Standard Model physics searches but also important for precise measurements of the top-quark fiducial cross sections and top-quark properties at the LHC.


Introduction
The top quark, discovered by the CDF and D0 experiments at Fermilab more than 20 years after its existence was postulated to explain the observed CP violations in kaon decays, is the heaviest elementary particle in the Standard Model (SM) of particle physics. Due to its large mass and its correspondingly short lifetime the top quark decays before hadronic bound states can be formed, thus, passing its spin information onto its decay products. With a mass of the order of the electroweak scale, the top quark Yukawa coupling to the Higgs boson is of the order of unity. This alone makes the top quark unique among the fermions and its potential to provide insights into physics beyond the SM (BSM) is anticipated. Various BSM models introduce modifications within the top quark sector, which can be constrained by precise measurements of the tt and tt + X cross sections, where X = H, j, γ, Z, W ± , tt. Examples include composite top quarks, Randall-Sundrum extra dimensions, models with coloured scalars or universal extra dimensions. Studies of top quark properties provide a unique environment for testing the SM and for hunting BSM physics. Investigations of the dynamics of the top quark pair production process in association with a hard photon, for example, directly probe the top quark electric charge and the structure of the ttγ coupling. Any deviation from the SM prediction of the measured observables could be an indication of BSM physics and might be linked to the production of an exotic (possibly heavier) top-like quark or the top quark with an anomalous electric dipole moment, see e.g. Refs. [1][2][3][4]. Good examples of such observables comprise the transverse momentum spectrum of the photon, p T,γ , and the azimuthal angle-rapidity distance between the photon and the softest b-jet, ∆R b 2 ,γ , [5,6].
First evidence for ttγ production has been reported by the CDF collaboration in pp collisions at the Tevatron with √ s = 1.96 TeV [7]. Observation was also announced by the ATLAS collaboration in pp collisions at the Large Hadron Collider (LHC) with √ s = 7 TeV [8]. Meanwhile, measurements have been carried out at the LHC by both ATLAS and CMS collaborations at √ s = 8 TeV [9,10]. For now, due to small available statistics, these measurements only comprise cross sections. However, with the second run of the LHC at √ s = 13 TeV and with increased luminosity more exclusive observables and various properties of the top quark can be scrutinised.
On the theory side, first next-to-leading order (NLO) QCD calculations for ttγ have been performed for on-shell top quarks [11][12][13]. Recently, even NLO electroweak corrections have been completed [14]. Computations in the approximation of stable tops give a general idea about the size of higher order effects. However, they can not provide a reliable description of top quark decay products or the magnitude of NLO corrections when specific cuts are imposed on the final states. For a realistic analysis, not only higher order effects to ttγ production need to be included but also radiative decays of top quarks have to be incorporated. This has been (partially) achieved by means of parton showers through matching fixed order NLO QCD predictions for ttγ with parton shower programs via the Powheg method [15,16], albeit omitting photon emissions in the parton shower evolution and tt spin correlations [17]. A more sophisticated approach has been employed in [18], where NLO QCD corrections to production and decays in the so called narrow width approximation (NWA) have been calculated. Non-factorisable QCD contributions, however, that imply a cross talk between production and decays of top quarks and which require going beyond the NWA, have been so far neglected. Such contributions are formally suppressed, i.e. O(Γ t /m t ) ≈ 0.8%, where Γ t , m t are the top quark width and mass respectively. They proved to be small in the inclusive cross section. Nonetheless, they can be strongly enhanced in case of exclusive observables that are crucial for new physics searches. The lack of any evidence of BSM at the LHC has put known new physics scenarios under significant strain. Our attention is focused now on precision physics and indirect searches aiming at deviations from SM predictions in precision observables. To probe more subtle BSM effects also in ttγ production, state of the art theoretical predictions for this process are of vital importance.
In this paper, we calculate for the first time the NLO QCD corrections to the fully realistic final state pp → e + ν e µ −ν µ bbγ. We consistently take into account resonant and non-resonant top quark and W gauge boson contributions and interference effects among them. Our theoretical predictions are presented in the form of the fully flexible Monte Carlo program. Thus, various observables and cuts can be explored and their usefulness can be demonstrated in realistic Monte Carlo simulations. The final results are provided as the Ntuple files [19]. Specifically, they are stored in the form of the modified Les Houches event files [20] and ROOT files [21] that might be directly employed in experimental studies at the LHC.
The article is organised as follows. In Section 2 we describe the details of our calculation. Input parameters and cuts to simulate detector response are summarised in Section 3. Numerical results for the integrated and differential cross sections for the LHC Run II energy of 13 TeV for two renormalisation (µ R ) and factorisation (µ F ) scale choices are presented in Sections 4 and 5, respectively. The theoretical uncertainties of the total cross sections and various differential cross sections, that are associated with neglected higher order terms in the perturbative expansion and with different parametrisations of the parton distribution functions (PDFs), are also given there. Finally, in Section 6 our conclusions are given.

Computational Framework
At leading order (LO) in the perturbative expansion, the pp → e + ν e µ −ν µ bbγ final state is produced via the scattering of either two gluons or a quark and the corresponding anti-quark gg → e + ν e µ −ν µ bbγ , qq → e + ν e µ −ν µ bbγ , where q stands for u, d, c, s. In total, the gg → e + ν e µ −ν µ bbγ subprocess comprises 628 Feynman diagrams and the qq → e + ν e µ −ν µ bbγ subprocess has only 346 tree level Feynman diagrams. Even though we do not use Feynman diagrams to obtain matrix elements for each subprocess, we provide such information here in order to shed some lights on the complexity of the process at hand. A few examples of Feynman diagrams contributing at O(α 2 s α 5 ) for the gg initial state are presented in Fig. 1. The calculation of scattering amplitudes is performed by means of an automatic offshell iterative algorithm [33], which is implemented within the Helac-Dipoles package [34] and the Helac-Phegas Monte Carlo (MC) program [35]. The latter framework has been used to cross check our LO results. For the phase-space integration depending on the MC framework Phegas [36], Parni [37] and Kaleu [38] have been employed.
At NLO, virtual corrections are obtained from the interference of the one-loop diagrams with the tree level amplitude. They might be classified into self-energy, vertex-, box-, pentagon-, hexagon-and heptagon-type corrections. For the gg dominant production channel we have 36032 Feynman diagrams at one-loop, among these 90 are heptagons and 958 are hexagons. The latter numbers have been obtained with the help of Qgraf [39]. Virtual corrections are evaluated in d = 4 − 2 dimensions in the 't Hooft-Veltman version of the dimensional regularisation and using the Feynman gauge for gauge bosons. The singularities coming from infrared divergent pieces are canceled by the corresponding ones arising from the counter-terms of the adopted subtraction scheme integrated over the phase space of the unresolved parton. The finite contributions of the loop diagrams are evaluated numerically in d = 4 dimensions. To ensure numerical stability of our calculations we have used the Ward identity test. On-shell transversality of gluon amplitudes has been checked up to the one loop level for every phase space point. About 6% of events, that fail the gauge-invariance check, have been recomputed in quadruple precision. For qq → e + ν e µ −ν µ bbγ partonic subprocess at O(α 3 s α 5 ) there are no gluons as external particles. Since unstable electroweak bosons are treated in the fixed-width scheme, the photon Ward identity test could not be applied straightforwardly. Instead the scale test, as introduced in Ref. [40], has been performed. It is based on the simple observation that the momenta can be rescaled and the amplitude can be recalculated and compared to the original one. In this case higher precision has also been used to recompute 0.15% of events that did not pass the test. Another cross check that we have performed comprises a verification of the cancelation of infrared poles. We compute the virtual corrections using Helac-1Loop [41] and CutTools [42], which are both parts of the Helac-NLO Monte Carlo framework [43]. The CutTools code contains an implementation of the OPP method for the reduction of one-loop amplitudes at the integrand level [44][45][46]. For unstable top quarks the complex mass scheme, as described in Refs. [47,48], is used. At the one loop level the appearance of Γ t = 0 in the propagator requires the evaluation of scalar integrals with complex masses, which is supported by the OneLOop program [49] employed in our studies. For consistency, mass renormalisation for the top quark is also done by applying the complex mass scheme in the well known on-shell scheme. Re-weighting techniques, helicity and colour sampling methods are employed in order to optimise the performance of our system.
The real emission corrections to the LO process arise from tree-level amplitudes with one additional parton, i.e. an additional gluon, or a quark anti-quark pair replacing a gluon. All possible subprocesses contributing to the real emission part are shown in Table  1. The number of Feynman diagrams corresponding to the subprocesses under scrutiny is also given to underline the complexity of the calculations. The following three subprocesses qg → e + ν e µ −ν µ bbγq,qg → e + ν e µ −ν µ bbγq and qq → e + ν e µ −ν µ bbγg are related by crossing symmetry. The singularities from soft and/or collinear parton emissions are isolated via subtraction methods for NLO QCD calculations: the commonly used Catani-Seymour dipole subtraction [34,50,51], and a fairly new Nagy-Soper subtraction scheme [52], both implemented in the Helac-Dipoles software. Specifically, Helac-Dipoles implements the massless dipole formalism of Catani and Seymour, as well as its massive version for arbitrary helicity eigenstates and colour configurations of the external partons. The Nagy-Soper subtraction scheme makes use of random polarisation and colour sampling of the external partons. An overall performance of this scheme has been assessed in Ref. [52] where a detailed comparison with results based on the Catani-Seymour dipole subtraction scheme has been carried out for a collection of processes. Thus, in Table 1 we only compare the total number of subtraction terms that need to be evaluated in both schemes. In each case, three times fewer terms are needed in the Nagy-Soper subtraction scheme compared qg → e + ν e µ −ν µ bbγq 2344 15 5 qg → e + ν e µ −ν µ bbγq 2344 15 5 qq → e + ν e µ −ν µ bbγg 2344 15 5 Table 1. The list of partonic subprocesses contributing to the subtracted real emission for the pp → e + ν e µ −ν µ bbγ + X process. Also shown are the number of Feynman diagrams, as well as the number of Catani-Seymour and Nagy-Soper subtraction terms.
to the Catani-Seymour scheme. The difference corresponds to the total number of possible spectators, which are relevant in the Catani-Seymour case, but not in the Nagy-Soper case. A phase space restriction (α max ) on the contribution of the subtraction terms, as proposed e.g. in Refs. [53][54][55][56][57][58], is included for both subtraction cases. We consider two choices, namely α max = 1, that corresponds to the original formulation of the Catani-Seymour and Nagy-Soper subtraction scheme, as well as α max = 0.01. In case of the Nagy-Soper subtraction scheme, which was our default scheme used for the calculations, we have checked that the final results for the sum of real radiation and integrated dipoles were independent of the α max choice. We have further cross checked that results for the real emission part are in agreement with results obtained with the Catani-Seymour dipole subtraction scheme. For the real correction part, we also adopt the Kaleu phase-space generator that is equipped with additional, special channels that proved to be important for phase-space optimisation.

Input Parameters and Cuts
In the following we present predictions for pp → e + ν e µ −ν µ bb γ + X production at O(α 3 s α 5 ) for the LHC Run II energy of √ s = 13 TeV. We consider decays of weak bosons to different lepton generations only, thus, neglecting the interference effects. However, the difference between the LO cross sections for pp → e + ν e µ −ν µ bb γ + X and pp → e + ν e e −ν e bb γ + X is at the per mille level, thus, the simplification is very well motivated. The complete cross section for pp → + ν −ν bbγ, with ± = e ± , µ ± , can be obtained by multiplying results presented in the following by a factor of 4. The Cabibbo-Kobayashi-Maskawa mixing of the quark generations is neglected and the following SM parameters are used The top quark width is calculated according to [59] at the scale m t . All other quarks, including b quarks as well as leptons, are assumed to be massless. Leptonic W gauge boson decays do not receive NLO QCD corrections. To include some effects of higher order corrections for the gauge boson widths, the NLO QCD values of the corresponding W width are used for LO and NLO matrix elements. The electroweak coupling is derived from the Fermi For our setup we have α Gµ ≈ 1/132. In the G µ -scheme electroweak corrections related to the running of α Gµ and to the ρ parameter (proportional to m 2 t /m 2 W ), are incorporated. By parametrising the lowest order in terms of G µ a large part of these universal electroweak corrections is absorbed. To describe the emission of the hard (real) photon, however, we use the α(0)−scheme with α ≡ α(0) = 1/137. Consequently, the prediction for the ttγ cross section is decreased by more than 3%. Based on the fact that relative NLO EW corrections to the on-shell ttγ production at the 13 TeV LHC are negative and of the order of 2% [14] we believe that this is a more consistent approach compared to employing α Gµ . In the first step we use kinematic-independent scales µ R = µ F = µ 0 with the central value µ 0 = m t /2 rather than µ 0 = m t . Even though the mass of the heaviest particle appearing in the process seems to be a more natural option, this scale choice is motivated by the fact that tt production at the LHC is dominated by t-channel gluon fusion, which favours smaller values of the scale. Additionally, the contributions beyond NLO that include the resummation of next-to-leading logarithmic soft gluon effects (NLO+NLL) are smaller for µ 0 = m t /2 than for µ 0 = m t , as we have explicitly checked with the help of the Top++ program [60]. Taking into account that photon emission is not a QCD effect this picture should not change for pp → ttγ production. With the goal of stabilising shapes in the high p T regions, that are relevant for the new physics searches, we have explored a dynamical choice for µ R and µ F . Kinematic-dependent scales should help to achieve flatter differential K-factors, thus, to describe more appropriately regions of the phase-space far away from the tt threshold. For the process at hand, we explored several possibilities and decided in the end to consider the following dynamical scale µ R = µ F = µ 0 = H T /4 where H T is the total transverse momentum of the system, which we have defined as where b 1 and b 2 are bottom-jets (not bottom quarks) and p miss T is the total missing transverse momentum from escaping neutrinos. The theoretical uncertainty is estimated with independent scale variation µ R = µ F , subject to the additional restriction 0 In practise such a restriction amounts to consider the following pairs Consequently, the minimum and maximum of the resulting cross section is chosen. Let us mention here that while calculating the scale dependence for the NLO cross section we keep Γ NLO t fixed independently of the scale choice. For two scales µ = µ 0 /2 and µ = 2µ 0 with µ 0 = m t /2 the change in the value of Γ NLO t is smaller than ±1.2%. The error introduced by this treatment is, however, of higher orders. We have checked that for the simpler case, i.e. for the pp → e + ν e µ −ν µ bb + X process, the variation in Γ NLO t has introduced deviations in the total cross section up to ±1.5% only [23]. Let us further note that in Ref. [24] a similar procedure has been discussed. In this paper the mismatch between the scale used in partial and total top quark decay widths has been compensated by the so-called partial width correction. The latter has been studied for total cross section for pp/pp → e + ν e µ −ν bb + X within cuts for Tevatron and LHC at different center of mass system energies both with fixed and dynamical scales. In the case of the LHC at √ s = 14 TeV, for example, these partial width corrections amounted to 1% − 3% depending on the scale choice employed. At NLO (LO) in QCD we employ CT14nlo (CT14llo) [61] PDFs and describe the running of the strong coupling constant α s with two-loop (one-loop) accuracy. The calculations are performed in the so-called 5 flavour scheme, however, the b-initiated contributions are not taken into account due to their numerical insignificance. To be more precise already at LO they are below 0.1%. All final-state partons with pseudorapidity |η| < 5 are recombined into infrared-safe jets via the anti−k T jet algorithm [62]. The cone size and jet resolution parameter R is set to R = 0.4. We require exactly two b-jets, one photon, two charged leptons and missing transverse momentum, p miss T . The hard photon is defined with p T, γ > 25 GeV and |y γ | < 2.5. To avoid QED collinear singularities in photon emission, caused by q → qγ splittings, a separation between quark and photon is required. Since distinguishing between quark and gluon jets is impossible on the experimental side, at the same time a separation between photons and gluons is induced as well. As a consequence, at a given photon p T an angular restriction on the soft gluon emission phase-space is introduced. Thus, soft divergences in the real emission part are different from those in the virtual correction impairing the cancelation of infrared divergences. To ensure soft and collinear safety we use a modified cone approach as described in Ref. [63], which implements a (smooth) isolation condition treating quarks and gluons the same way. With the isolation cone of R γj = 0.4 for each parton i we evaluate ∆R γi between this parton and the photon. We reject the event unless the following condition is fulfilled for all R ≤ R γj , where E T, i is the transverse energy of the parton i and E T, γ is the transverse energy of the photon. Jets reconstructed inside the cone size R γj are not subject to additional selection criteria. We apply the following inclusive cuts to simulate detector response where stands for µ − , e + . We set no restriction on the kinematics of the extra jet.  4 Results for the LHC Run II energy of 13 TeV for the fixed scale choice

Integrated cross section and its scale dependence
With the input parameters and cuts specified above, we arrive at the following predictions  At the central scale µ 0 = m t /2, the gg channel dominates the total LO pp cross section by 79%, followed by the qq +qq channel with 21%. Photons are, therefore, predominantly radiated off the top quark and top quark decay products. The full pp cross section receives negative and moderate NLO corrections of 10%. The theoretical uncertainties resulting from scale variations, where µ R and µ F have been varied independently, and taken in a very conservative way as a maximum of the lower and upper bounds are 35% at LO and 14% at NLO. Thus, a reduction of the theoretical error by a factor of 2.5 is observed. Should we instead vary µ R and µ F simultaneously, up and down by a factor of 2 around µ 0 , the uncertainties would remain unchanged. This is due to the fact that the scale variation is driven solely by the changes in µ R . In the case of truly asymmetric uncertainties, however, it is always more appropriate to symmetrise the errors. After symmetrisation the scale uncertainty at LO is assessed to be instead of the order of 30%. After inclusion of the NLO QCD corrections, the scale uncertainty is reduced down to 7%. The graphical presentation of the behaviour of LO and NLO cross sections upon varying the scale by a factor ξ ∈ {0.125, . . . , 8} is shown in Fig. 2. At LO the individual contributions of the partonic subprocesses are additionally presented. The final scale dependence of the NLO cross section as emerged out of the two contributions (the virtual plus the LO part and the real emission part) is also depicted in Fig. 2. Of course, the separation is entirely unphysical, but well defined once we state that we use the 't Hooft-Veltman version of the dimensional regularisation, with the integrals as defined in the OneLOop library. Next, we have checked the dependence on the parameters introduced in the photon isolation procedure. Specifically, the general photon isolation formula is given by with two additional parameters γ and n. The default choice, which should guarantee moderate corrections, is γ = 1 and n = 1, see Ref. [63]. Nevertheless, both γ and n a priori can have arbitrary values. We have recalculated the subtracted real emission part of the NLO results with a different choice, namely γ = 1/2 and n = 1/2. Within the integration errors our new results have agreed with the old ones. Thus, NLO QCD results for the pp → e + ν e µ −ν µ bbγ + X production process are not sensitive to moderate changes in values of γ and n.
In the subsequent step, the size of the top quark non-factorisable corrections has been estimated for the total cross section. To achieve this the full result has been compared with the result in the NWA. The latter has been obtained by rescaling the coupling of the top quark to the W boson and the b quark by several large factors, as described in Ref. [23], to mimic the limit Γ t → 0 when the scattering cross section factorizes into on-shell production and decay. The top quark non-factorisable corrections for the pp → e + ν e µ −ν µ bbγ + X production process amount to 1.5% (2.5%) for LO (NLO). They are consistent with the expected uncertainty of the NWA, which is of the order of O(Γ t /m t ).
Coming back to the theoretical uncertainties, we note that, another source of theoretical uncertainties comes from the PDF parametrisation. To that end, we have recomputed NLO QCD corrections to the pp → e + ν e µ −ν µ bbγ +X production process with different PDF sets. Following recommendations of PDF4LHC for the usage of PDFs suitable for applications at the LHC Run II [64] we employ additionally to the CT14 PDF set the MMHT14 PDF set [65] and NNPDF3.0 [66]. Let us say here, that differences coming from NLO results  Figure 3. Differential distributions as a function of the transverse momentum of the hard photon, p T,γ and the rapdity-azimuthal angle separation between the photon and the softer b-jet, ∆R b2,γ , for µ F = µ R = µ 0 = m t /2. The LO and the NLO CT14 PDF sets are employed. The upper panels show absolute LO and NLO predictions together with corresponding uncertainty bands. The lower panels display the differential K-factor together with the uncertainty band and the relative scale uncertainties of the LO cross section.
for various PDF sets are comparable (usually even higher) to the individual estimates of PDF systematics. We have checked that this is the case for the similar process, namely for pp → e + ν e µ −ν µ bbj + X production [30]. In this paper, we take the PDF uncertainties to be the difference between our default PDF set (CT14) and the other two PDF sets considered (MMHT14 and NNPDF3.0). Our findings for MMHT14 and NNPDF3.0 PDF can be summarised as follows σ NLO pp→e + νeµ −ν µbbγ (MMHT14, µ 0 = m t /2) = 7.49 fb , σ NLO pp→e + νeµ −ν µbbγ (NNPDF3.0, µ 0 = m t /2) = 7.72 fb .

(4.3)
The PDF uncertainties for the process under scrutiny are, therefore, given by +0.05 fb (1%) for the MMHT14 PDF set and +0.28 fb (4%) for NNPDF3.0. Our result for the integrated cross section at NLO in QCD with the CT14 PDF set and for µ 0 = m t /2 is given by Taken in a very conservative way, the PDF uncertainties are of the order of 4% (to be compared to the theoretical uncertainties of 14% from the scale dependance). After symmetrisation they are reduced down to 2% (to be compared to 7%). Overall, the PDF uncertainties are well below the theoretical uncertainties due to the scale dependence. The latter remain the dominant source of the theoretical systematics.  Figure 4. Differential distributions as a function of the separation in the rapidity-azimuthal angle plane between the hard photon and the hardest b-jet, ∆R b1,γ , as well as the separation between the hard photon and the hardest and softer lepton, ∆R 1,γ and ∆R 2,γ for µ F = µ R = µ 0 = m t /2. The LO and the NLO CT14 PDF sets are employed. The upper panels show absolute LO and NLO predictions together with corresponding uncertainty bands. The lower panels display the differential K-factor together with the uncertainty band and the relative scale uncertainties of the LO cross section.

Differential cross sections
While the size of higher order corrections to the total cross section is certainly interesting, it is crucial to study the corrections to differential distributions. In Fig. 3 we present representative differential distributions, that are relevant for BSM searches [5,6]. We display p T of the hard photon and ∆R b 2 ,γ between the hard photon and the softer b-jet. The upper panels show the distributions themselves and their scale dependence. The lower panels reveal the differential K-factor with its error and the relative scale uncertainties of the LO cross section. To be more precise we plot K NLO (µ) = (dσ NLO (µ)/dX)/(dσ LO (µ 0 )/dX)  Figure 5. Differential distributions as a function of p T of the hardest and the softer b-jet as well as the hardest and the softer lepton for µ F = µ R = µ 0 = m t /2. The LO and the NLO CT14 PDF sets are employed. The upper panels show absolute LO and NLO predictions together with corresponding uncertainty bands. The lower panels display the differential K-factor together with the uncertainty band and the relative scale uncertainties of the LO cross section. and K LO (µ) = (dσ LO (µ)/dX)/(dσ LO (µ 0 )/dX) where µ 0 = m t /2 is the central value of the scale and X denotes the observable that is scrutinised. Higher order corrections have strongly altered the shape of ∆R b 2 ,γ where corrections range from −29% to +122%, causing distortions of up to 150%. Similar effects have been noticed for other observables, most notably for other angular observables shown in Fig. 4. Among others the most affected by higher order corrections are the separation in the rapidity-azimuthal angle plane between the hard photon and the hardest b-jet, ∆R b 1 ,γ (NLO corrections from −24% to +93%), as well as the separation between the hard photon and the hardest or softer charged lepton, ∆R 1 ,γ (NLO corrections ranging from −25% to +91%) and ∆R 2 ,γ (NLO corrections ranging from −16% to +132%). In each case the large differential K-factor for ∆R 4 is associated with photon emission from initial state quark from the qg + gq partonic subprocess, where q stands for quark and antiquark. Such a contribution appears only starting at NLO and adds significantly at large ∆R. Let us mention here, that the qg + gq channel contribution to σ NLO pp→e + νeµ −ν µbbγ is estimated at the level of 29%. Moreover, due to the leading order like nature of the contribution also the scale dependence in this region is enlarged. Let us additionally note here, that emission of the photon from the charged lepton leads to collinear enhancement at small values of the separation between the photon and the softer charged lepton, ∆R 2 ,γ as can be clearly observed in Fig. 4. Moreover, in the case of the separation between the photon and the softer b-jet, R b 2 ,γ , depicted in Fig. 3, events are produced over a wide range of ∆R b 2 ,γ values rather than in the back-to-back configuration. This confirms the findings of Ref. [18] that photon radiation off top quark decay products yields a significant contribution to the cross-section.
In case of p T,γ the differential K-factor is rather constant only in the plotted range from −8% to −18%. Thus, p T of the photon is more stable against higher order corrections and hence better suited for BSM searches. Nevertheless, both observables require higher order calculations to be properly described. In view of ongoing indirect searches for BSM physics, where the emphasis is looking for small deviations from the most accurate SM predictions, such state of the art results are indispensable.
Finally, in Fig. 5 we present dimensionful observables. Specifically, we display the transverse momentum distributions of the hardest and the softer b-jet and charged lepton. In all cases negative and large higher order QCD corrections can be detected. In the high p T regions they amount to −38%, −53%, −43% and −76% for p T, b 1 , p T, b 2 , p T, 1 and p T, 2 respectively. Moreover, the NLO error bands do not fit within the LO ones as one would expect from a well-behaved perturbative expansion. Thus, the fixed scale choice does not ensure a stable shape when going from LO to NLO for these observables. Through the implementation of a dynamical scale, the large discrepancies between the shapes of these distributions at NLO and LO should disappear. Thus, in the next step we shall examine NLO results for µ R = µ F = µ 0 = H T /4 with the goal of stabilising differential K-factor, i.e. decreasing NLO QCD corrections in the tails, for p T, b 1 , p T, b 2 , p T, 1 and p T, 2 while keeping the behaviour of K almost unchanged for p T, γ .

Theoretical uncertainties for differential cross sections
At this point we would like to estimate theoretical uncertainties inherent in our LO and NLO differential cross sections as obtained with µ 0 = m t /2 and the CT14 PDF set. The scale uncertainties are again estimated conservatively by scanning bin by bin values of the lower and upper bounds and by choosing the maximal number. To get a general idea about the size of theoretical errors we quote here only this maximal value. In this way, for the transverse momentum distribution of the hard photon we have obtained theoretical errors up to ±40% at LO and up to ±22% at NLO. For dimensionful observables these maximal values come from the high p T regions. A similar pattern can be seen for the angular separation between the hard photon and b-jet or the charged lepton. Specifically, for ∆R b 2 ,γ we have ±40% at LO to be compared with ±33% at NLO and for ∆R b 1 ,γ is ±42% at LO versus ±28% at NLO. In case of the charged lepton the situation is very similar as we have estimated the theoretical error at the level of ±42% (±40%) at LO and ±28% (±27%) at NLO for ∆R 1 ,γ (∆R 2 ,γ ). Thus, in all above mentioned cases a reduction by a factor of 1.5 − 2 is achieved by increasing the order in perturbative expansion. However, in case of transverse momentum distributions of b-jets and leptons the picture has changed and there is a large residual scale dependence in these observables even at NLO. Actually, for the p T,b 1 distribution the theoretical error is at the same level independently of the perturbative order and amounts to ±46%. For the p T, 1 and p T,b 2 at NLO the theoretical error is larger than at LO, respectively ±56% and ±78%. Finally, for the last plotted observable, i.e. p T, 2 huge uncertainties of the order of ±186% can be seen. This clearly tell us that µ 0 = m t /2 is not equipped to properly describe tails of p T distributions even at NLO. Many of these features can be improved by performing NLO computation with the kinematic-dependent choices of the scales.
Lastly, we have examined PDF uncertainties for the differential cross sections with the fixed scale choice. For all observables that we have studied PDF uncertainties are negligible in comparison to the theoretical uncertainties from the scale dependence. As an example we show in Fig. 6 NLO differential distributions as a function of the transverse momentum of the hard photon and the azimuthal angle-rapidity distance between the hard photon and the softest b-jet. The upper panels present the NLO predictions for three different PDF sets at the central scale value µ R = µ F = µ 0 = m t /2. In addition to the CT14 PDF set, we employ the MMHT14 and NNPDF3.0 PDF sets. The lower panels of Fig. 6 give the ratio of the MMHT14 (NNPDF3.0) PDF set to CT14.
To summarise this part, for pp → e + ν e µ −ν µ bbγ + X production at the LHC Run II with √ s = 13 TeV with our selection of cuts and input parameters, the PDF uncertainties are insignificant both at the level of total and differential cross sections once contrasted with theoretical errors from the scale dependence. Let us note at this point, however, that additional theoretical effects should be examined for the process at hand. These include among others NLO electroweak effects, the size of which has to be estimated and compared to the size of NLO QCD effects. Moreover, dedicated analyses of complete NLO QCD off-shell effects of the top quark at the differential level have to be carried out. We leave both aspects for future studies. Thus, from now on we shall concentrate only on theoretical uncertainties from the scale dependence.  Fig. 7. Moreover, the variation of µ R (µ F ) with fixed µ F (µ R ) is presented in Fig. 7 as well. Here qualitatively our findings remain the same as for the fixed scale choice.  We start with differential distribution for the transverse momentum of the hard photon that is displayed in Fig. 8. For the dynamical scale choice of µ 0 = H T /4 positive corrections up to 13% are obtained. We can also notice that the NLO error band as calculated through the scale variation is within the LO error band as it should be for a well behaved observable that is described by the perturbative expansion in α s . For the dimensionless observable ∆R b 2 ,γ , that is also shown in Fig. 8, the size of NLO corrections has been moderately reduced. The higher order corrections range now from −12% up to +116%. Thus, the shape distortion up to 128% has been obtained for this observable, which should be compared with 150% for the fixed scale choice. Other dimensionless observable are presented in Fig. 9. For the differential distribution as a function of the separation in the rapidity-azimuthal angle plane between the hard photon and the b-jet or lepton, i.e. ∆R b 1 ,γ , ∆R 1 ,γ and ∆R 2 ,γ we have acquired NLO QCD corrections in the following range {−11%, +90%}, {−10%, +90%} and {−5%, +130%} respectively. In each case shape distortions have been decreased by about 15% − 20%.
Finally, we have reexamined dimensionful observables. Specifically, transverse momentum distributions of the hardest and the softer lepton as well as transverse momentum distributions of the hardest and the softer b-jet. They are given in Fig. 10. Higher order corrections in high p T regions have been substantially reduced. For the transverse momentum distribution of the hardest b-jet we have attained +19% instead of −38% and for the softer b-jet −16% to be compared with −53% for the fixed scale choice. The same pattern can be noticed for the p T differential cross section of the hardest (the softer) lepton. In details, for the hardest one we have obtained +8% as a substitute to −43% whereas for the softer charged lepton −30% rather than −76%.
To summarise this part, the validity of the proposed dynamical scale µ 0 = H T /4, that is blind to the underlining top quark resonance history, is confirmed. The size of NLO QCD corrections to all presented observables has been reduced. Moreover, this judicious choice of the scale has allowed us to obtain nearly constant K-factors in all dimensionful distributions that we have studied.

Theoretical uncertainties for differential cross sections
As a final step we have examined theoretical uncertainties for differential cross sections for the dynamical scale choice µ 0 = H T /4. The µ 0 = m t /2 scale choice has proved to be inadequate for the modelling of various differential distributions and more importantly for the estimation of their theoretical errors in the high p T regions. The latter phase space regions are simply not very sensitive to the threshold contributions for the ttγ production that are well described by the fixed scale choice. For each considered observable we have  Figure 9. Differential distributions as a function of the separation in the rapidity-azimuthal angle plane between the hard photon and the hardest b-jet, ∆R b1,γ , as well as the separation between the hard photon and the hardest and softer charged lepton, ∆R 1,γ and ∆R 2,γ for µ F = µ R = µ 0 = H T /4. The LO and the NLO CT14 PDF sets are employed. The upper panels show absolute LO and NLO predictions together with corresponding uncertainty bands. The lower panels display the differential K-factor together with the uncertainty band and the relative scale uncertainties of the LO cross section.
observed reduced theoretical errors as compared to the µ 0 = m t /2 scale choice. The effect is more pronounce in the case of dimensionful observables in the high p T regions. Thus, for example we can see from Fig. 8 where the differential cross section as a function of the transverse momentum of the hard photon is plotted, that theoretical error is now up to ±8% at NLO (up to ±36% at LO) as compared to ±22% at NLO (±40% at LO) for µ 0 = m t /2. When considering p T distribution of the hardest and the softer b-jet, depicted in Fig. 10, the theoretical error at NLO is reduced from ±47% and ±78% down to ±10% and ±18% respectively. The most dramatic effect can be seen in the case of p T distribution of the hardest and the softer lepton, also given in Fig. 10. In that case instead of theoretical  Figure 10. Differential distributions as a function of p T of the hardest and the softer lepton as well as the hardest and the softer b-jet for µ F = µ R = µ 0 = H T /4. The LO and the NLO CT14 PDF sets are employed. The upper panels show absolute LO and NLO predictions together with corresponding uncertainty bands. The lower panels display the differential K-factor together with the uncertainty band and the relative scale uncertainties of the LO cross section. errors up to ±56% and ±186% we have received theoretical errors up to only ±7% and ±31% respectively.
To recapitulate this part, the dynamical scale µ 0 = H T /4, which has been considered in our analysis, has proven to be very effective in stabilising the perturbative convergence in the phase space regions far away from the 2m t threshold and in providing small theoretical uncertainties as estimated by the scale variation. For all considered observables the latter are below 10% − 30%.

Summary and Outlook
We have presented the first complete higher order predictions for the pp → ttγ process in the di-lepton channel for the LHC run II energy of √ s = 13 TeV. With our inclusive cuts and for µ R = µ F = µ 0 = m t /2, NLO predictions reduced the unphysical scale dependence by a factor of 2.5 and lowered the total rate by about 10% compared to LO predictions. The theoretical uncertainty of the NLO cross section as obtained from the scale dependence has been estimated to be 14%. By comparison the PDF uncertainties are negligible at the level of 4% only. On the other hand, for µ R = µ F = µ 0 = H T /4 the full pp cross section has received positive and small NLO QCD corrections of 2.5%. Additionally, the inclusion of higher order effects has reduced the theoretical error by a factor of 5.5. Specifically, the theoretical uncertainties due to scale dependence are now at the 6% level only, however, they are still larger than the PDF uncertainties.
Even though NLO QCD corrections to the total cross section vary from moderate to small depending on the scale choice their impact on differential distributions is much larger. Independently of the scale choice for some dimensionless observables shape distortions of more than 100% have been observed. The prominent example comprises the differential cross section as a function of the separation in the rapidity-azimuthal angle plane between the hard photon and the softer b-jet, ∆R b 2 ,γ . In the case of this observable, which is relevant for new physics searches, shape distortions up to 150% (128%) have been obtained for µ 0 = m t /2 (µ 0 = H T /4). For the dimensionful observables presented in this paper, however, the dynamical scale choice has helped to obtain almost flat differential K-factors as well as to stabilise the high p T regions, which are very poorly described by NLO results with the fixed scale choice. Also in the case of differential observables the PDF uncertainties have been examined. Similarly to the total cross section case their size is negligible when comparing to scale uncertainties. We repeat at this point that additional theoretical effects should be investigated. Among others the size of NLO electroweak effects has to be calculated for the pp → e + ν e µ −ν µ bbγ cross section and for various differential cross sections. We plan to include such effects in a future publication.
In addition, the size of the top quark off-shell effects for the total cross section has been estimated to be 2.5%. Their influence on differential distributions and extraction of the SM parameters, however, might be much stronger, as has already been shown in case of tt and ttj production [67][68][69]. Again, we leave such studies for the future.
Our theoretical predictions are stored in the form of the Ntuples files and are available upon request. Specifically, they are stored in the form of modified Les Houches event files and ROOT files, that might be directly employed in experimental studies at the LHC. They can be used for example to change kinematical cuts or to define new observables. The latter can be obtained without need of any additional rerunning of the code. Moreover, any change in the renormalisation or factorisation scale choice or in the PDF set can be accommodated by simple reweighting of these files. Thus, they can be employed to study broad phenomenological aspects of top quark physics at the LHC.