Constraining anomalous Higgs boson couplings to virtual photons

We present a study of Higgs boson production in vector boson fusion and in association with a vector boson and its decay to two vector bosons, with a focus on the treatment of virtual loops and virtual photons. Our analysis is performed with the JHU generator framework. Comparisons are made to several other frameworks, and the results are expressed in terms of an effective field theory. New features of this study include a proposal on how to handle singularities involving Higgs boson decays to light fermions via photons, calculation of the partial Higgs boson width in the presence of anomalous couplings to photons, a comparison of the next-to-leading-order electroweak corrections to effects from effective couplings, and phenomenological observations regarding the special role of intermediate photons in analysis of LHC data in the effective field theory framework. Some of these features are illustrated with projections for experimental measurements with the full LHC and HL-LHC datasets.


I. INTRODUCTION
The large amount of data analyzed by the ATLAS and CMS experiments on the Large Hadron Collider (LHC) is consistent with the predictions of the standard model (SM) of particle physics. Among the measurements performed, the discovery and characterization of the Higgs (H) boson have been crucial in completing the SM [1][2][3]. Yet, open questions remain, such as the low value of the H boson's mass, its Yukawa coupling hierarchy, the source of CP violation required for matter abundance, and the connection of the SM to other cosmological observations. Studies of electroweak production (VBF and V H) and decay (H → V V ) of the H boson probe HV V interactions over a large range of momentum transfer, which can expose possible new particles that couple through loops. Such electroweak processes lead to rich information in kinematic distributions of the H boson decay products and associated particles. Analysis of such distributions can shed light on the nature of the HV V interactions and has been discussed extensively . Such studies can be naturally performed within the effective field theory (EFT) framework, with examples of application to LHC data documented in Refs. [3,[46][47][48][49][50][51][52][53][54][55][56][57][58][59][60].
In our earlier studies using the JHU generator framework [19,23,28,38,45,61], we relied on dedicated Monte Carlo simulation, and demonstrated optimal discrimination, reweighting techniques, and analysis of a bosonic resonance with the most general anomalous couplings. We build upon this framework of the JHU generator and MELA analysis package with the goal of demonstrating its application to the H boson's interactions in electroweak processes with massless vector bosons, such as in HZγ, Hγγ, and Hgg vertices. Such couplings are generated in the SM through loops of SM particles. They also lead to divergence of fixed-order calculations for virtual γ * states when the fourmomentum squared q 2 γ * approaches zero. In the perturbative expansion, such terms are poorly defined at low values of q 2 . Some prior discussion of this effect can be found in Refs. [40][41][42][43][44]62].
We review the parameterization of anomalous H boson couplings in Section II and discuss applications of the JHU generator framework to EFT studies in Section III. In Section IV, the partial H boson width and production cross sections are calculated in the presence of anomalous couplings to massless vector bosons. In Section V, the treatment of the next-to-leading-order (NLO) electroweak (EW) effects is discussed. In Section VI, we make a proposal on how to handle singularities involving intermediate photons in the H boson decays. In Section VII, several phenomenological observations are made in application to LHC data.

II. PARAMETERIZATION OF ANOMALOUS INTERACTIONS
We start with the HV V scattering amplitude of a spin-zero boson H and two vector bosons V V with polarization vectors and momenta ε µ 1 , q µ 1 and ε µ 2 , q µ 2 . The amplitude is parameterized by where v is the vacuum expectation value, under the conventions ε 0123 = +1 and (q µ ) = (E, q). This amplitude represents the three possible tensor structures of the H boson's interaction with two vector bosons, with expansion of the terms up to q 2 i . By symmetry we have κ ZZ 1 = κ ZZ 2 , but we do not enforce κ W W for W ± bosons. Note that κ γγ 1 = κ γγ 2 = κ gg 1 = κ gg 2 = κ Zγ 1 = 0, while κ Zγ 2 /(Λ Zγ 1 ) 2 may contribute. The coupling κ V V 3 /(Λ V V Q ) 2 allows for scenarios which violate the gauge symmetries of the SM.
An effective HV V interaction may be generated by loops of fermions, in which case the couplings κ f andκ f describe the H boson interactions as whereψ f and ψ f are the Dirac spinors and m f is the fermion mass. For the SM fermions, (κ f ,κ f ) = (1, 0). The equivalent Lagrangian for H boson interactions with gauge bosons (in the mass eigenstate parameterization) reads A µν A µν +c γγ e 2 4 A µνÃ µν + c gg g 2 s 4 G a µν G aµν +c gg in accordance with Eq. (II.2.20) in Ref. [39], where e 2 = 4πα and g 2 s = 4πα s are the squared electromagnetic and strong coupling constants, respectively, and s w = 1 − c 2 w is the sine of the weak mixing angle. The covariant derivative used to derive this expression is D µ = ∂ µ − i e 2sw σ i W i µ − i e 2cw B µ [39,63]. We note that the convention ε 0123 = +1 defines the relative sign of the CP -oddc i and CP -even c i couplings [45], while the relative sign in front of the W i µ and B µ terms in the covariant derivative defines the sign of the Zγ couplings relative to the ZZ and γγ. The latter could be viewed as the sign of s w , if a different convention is adopted. 1 The generality of our amplitude parameterization allows us to uniquely represent each EFT coefficient in Eq. (3) by an anomalous coupling in Eq. (1).
Note that not every anomalous coupling in Eq. (1) has a corresponding term in the EFT Lagrangian of Eq. (3). For example, the term κ V V 3 /(Λ V V Q ) 2 is not gauge invariant and is not present in Eq. (3). Similarly, κ W W 1 = κ W W 2 due to charge symmetry.
So far we have discussed the H boson interactions without considering additional symmetries. The SU(3) × SU(2) × U(1) symmetry of the standard model effective field theory (SMEFT) [64][65][66][67] is a motivated framework which allows relating EFT operators. Not all of the EFT coefficients are independent when limiting the discussion to dimension-six interactions with this symmetry. The linear relations for the dependent coefficients can be found in Ref. [39] and they translate into relations amongst our anomalous couplings as follows: κ Zγ The Lagrangian for H boson interactions with gauge bosons can be written in the Warsaw basis [68] which preserves the SU(3) × SU(2) × U(1) symmetry of SMEFT. The relationship between operators in the Warsaw basis and the mass-eigenstate basis is discussed in Section III.

III. THE JHU GENERATOR FRAMEWORK AND THE EFT BASES
The JHU generator framework (JHUGen) includes a Monte Carlo generator and matrix element techniques for optimal analysis of the data. It is built upon the earlier developed framework of the JHU generator and MELA analysis package [19,23,28,38,45,61] and extensively uses matrix elements provided by MCFM [69][70][71][72][73]. The SM processes in MCFM are extended to add the most general scalar and gauge couplings and possible additional states. This framework includes many options for production and decay of the H boson, which include the gluon fusion, vector boson fusion, and associated production with a vector boson (V H) in both on-shell H and off-shell H * production [45]. In the off-shell case, interference with background processes or a second resonance is included. The processes with direct sensitivity to fermion Hff couplings, such as ttH, bbH, tqH, tW H, or H → τ + τ − , are discussed in Refs. [38,61].
The JHUGen framework was adopted in Run-I analyses using LHC data [2, 3, 46-50, 52, 74-77] and employed in recent Run-II measurements of the HV V anomalous couplings from the first joint analysis of on-shell production and decay [53,58], from the first joint analysis of on-shell and off-shell H boson production [57], for the first measurement of the CP structure of the Yukawa interaction between the H boson and top quark [78], in the search for a second resonance in interference with the continuum background [79,80], and in EFT approach to the HV V , Hgg, and Hff interactions [60].

A. EFT basis considerations
The framework is based on the amplitude parameterization in Eqs. (1) and (2). In order to simplify translation between different coupling conventions and operator bases, including the Higgs and Warsaw bases, within the JHU generator framework, we provide the JHUGenLexicon program, which includes an interface to the generator and matrix element library and can also be used for standalone or other applications [45]. The relationship of the amplitude parameterization to the mass eigenstate basis of the EFT formulations in Eq. (3) is performed through the simple linear relationship in Eq. (4). The functionality of this program is similar to Rosetta [81], but it is limited in scope to application to the H boson interactions and provides additional options to introduce certain symmetries or constraints, as illustrated below.
We count five CP -even and three CP -odd independent electroweak HV V operators, as well as one CP -even and one CP -odd Hgg operators in the mass-eigenstate basis in Section II. The same number of independent H boson operators exists in the Warsaw basis. The relationship between the six CP -even operators is quoted explicitly in Eq. (14) of Ref. [81]. This relationship is direct, with the exception of the δv parameter defined in the Warsaw basis in Eq. (15) of Ref. [81]. One could remove an extra parameter from transformation with constraints from precision electroweak data. For example, we can set ∆M W = 0 in Eq. (5), because M W is measured precisely. This allows us to express δv through the other HV V operators in the Warsaw basis. The JHUGenLexicon program provides such an option and the following studies in this paper will be presented with such a constraint. With the above symmetries and constraints, including ∆M W = 0, the translation between the Warsaw basis and the independent amplitude coefficients is where Λ is the scale of new physics, which we set to Λ = 1 TeV as a convention, and δg ZZ 1 is the correction to the SM value of g ZZ 1 = 2. According to Eq. (5), δg W W 1 = δg ZZ 1 , and the other dependent amplitude coefficients can be derived from Eqs. (6)(7)(8)(9).
A numerical example of the relationship between the C HX = 1 contribution of a single operator in the Warsaw basis and the couplings in the mass-eigenstate amplitude in Eq. (1) is shown in Table I, which corresponds to the reverse of Eq. (10).
B. Application to the VBF, V H, and H → V V processes One of the new features in this paper, compared to the earlier work, is the study of the g Zγ 2 , g γγ 2 , g Zγ 4 , and g γγ 4 anomalous couplings in electroweak production of the H boson. Their effect in the H → 4 process was studied with LHC data [48] and with phenomenological tools [25,35]. In the following we re-examine the H → 4 decay and investigate the VBF and V H processes. In the case of V H production, we consider three final states Z(→ ff )H, γ * (→ ff )H, and γH, and both qq or gg production channels, as all are affected by the HV V couplings of our interest. While the gluon fusion process formally appears at higher order in QCD, the large gluon parton luminosity at the LHC makes this channel interesting to examine. In this study, we only examine the operators affecting the H boson interactions in Table I and study their effect on the HV V couplings. Other operators, such as HZf f contact terms for example, are included in the JHUGen framework, but they are not the primary interest in this study because their existence would become evident in resonance searches and in electroweak measurements, without the need for H boson production. Moreover, such contact terms are equivalent to the combination of the κ ZZ 1 and κ Zγ 2 couplings if flavor universality is assumed [45]. Not only HV V interactions may be affected by the above operators in the processes under study. For example, the C HW B operator also affects the Zf f couplings. However, these Zf f couplings should be well constrained in electroweak measurements. For this reason, should one of the considered operators affect the Zf f interactions, we assume some other operators not affecting the direct H boson interactions must also contribute to bring the Zf f couplings to the SM values.
Numerical results of the relative contributions of operators to the H → V V → 4 , VBF, qq or gg → V (→ + − )H, and γH processes are shown in Appendix A. The general observations from Tables IV-VI is that the relative importance of the g Zγ 2 , g γγ 2 , g Zγ 4 , and g γγ 4 couplings changes between the processes. Taking the example of the C HW B operator, these couplings lead to an overwhelming contribution in the H → 4 process. However, their contribution in the VBF and V H processes is not significant and is especially tiny in the case of the V H process. These features will affect our ability to use different processes to constraint anomalous couplings with photons. We note that the V H process with V → + − includes both ZH and γ * H production mechanisms, where γ * leads to low-q 2 contributions in the m invariant mass, which can be observed in kinematic distributions.
The gg → ZH process has been shown to have no contributions of the two anomalous HV V tensor structures appearing in Eq. (1) in the triangular loop diagram [45]. Therefore, only the SM-like tensor structure with the g 1 and κ ZZ 1 couplings contributes to this diagram, as shown in Table VII. The off-shell photon does not couple to the triangular fermion loop either [45], and, therefore, the κ Zγ 2 coupling does not contribute. The box diagram is sensitive to the fermion couplings of the H boson, which we do not vary in this study of anomalous HV V interactions. As the result, the gg → ZH process features a rather limited set of EFT operators and we will not study this process in more detail in this paper, leaving further details to Ref. [45].
The γH production process has been largely neglected in analysis of LHC data. However, this process was used in the search for the H boson with anomalous couplings in e + e − production prior to the H boson discovery [82] and proposed in application to CP -even EFT operator constraints at the LHC [83,84]. From Table VIII, it is evident that only the g Zγ 2 , g γγ 2 , g Zγ 4 , and g γγ 4 couplings contribute, and this channel does not receive tree-level SM contributions. Because the photon is on-shell, it does not receive contribution from κ Zγ 2 either. This process is generated by the dimension-6 operators squared in the EFT expansion in combination with the EW loops generated by the SM particles. As an approximation to the SM production cross section, we use the calculation with the g Zγ,SM 2 and g γγ,SM 2 values calculated in Section V. These point-like couplings reproduce the SM decay width of the processes H → Zγ and γγ, respectively. Due to the off-shell V = Z/γ * in the process qq → V → γH, these point-like couplings are not expected to reproduce the full EW loop calculation in the SM, but they are expected to provide a good estimate, which we use as σ γH SM in Table VIII. The γH process may be of particular interest in isolating the CP -odd couplings g Zγ 4 and g γγ 4 in combination with CP -even couplings g Zγ 2 and g γγ 2 , which is complementary to the H → Zγ and γγ decays.

C. Kinematic distributions with EFT effects
The kinematic effects in the H → V V , VBF, and V H processes can typically be described with five angular observables and two invariant masses, or q 2 i of the two vector bosons, as illustrated in Fig. 1 [19,28,45]. The distributions of two of these angles, θ * and Φ 1 , are random for a spin-zero H boson, but are less trivial for a higherspin resonance or non-resonant production. In the following, we disentangle the relative contributions of the ZZ, W W , Zγ, and γγ intermediate vector-boson states to simulation with a given operator in the Warsaw basis. Such a decomposition reveals interesting kinematic effects and also allows us to validate the tools used for simulation of EFT effects and match their conventions.
We use the JHUGen program to generate several models which allow us to visualize the relative contributions of the mass eigenstates of the vector bosons. We also model the inclusive kinematic distributions with the SMEFTsim program [85] using MadGraph5 simulation [86]. Once the sign conventions are matched, as discussed in Section II, we find good agreement. A similar comparison with SM couplings using the Prophecy4f [87] and HAWK [88] generators is shown in Section V. Throughout this paper and unless otherwise noted, the calculations are performed FIG. 1: Three kinematic topologies of the H boson production and decay [45]: vector boson fusion at LO in QCD and EW, with the MS-mass for the top quark m t = 162.7 GeV, the on-shell mass for the bottom quark m b = 4.18 GeV, QCD scale µ = M H /2, α s = 0.1188, α = 1/128, s 2 w = 0.23119, G F = 1.16639 × 10 −5 GeV −2 [89], and the NNPDF 3.0 parton distribution functions [90].
In the H → 4 and V H processes, we require m > 1 GeV. In the VBF process, we apply the selection requirements m jj > 300 GeV, p jet T > 1 GeV, |η jet | < 5, ∆η jj > 1, ∆R jj > 0.3, q 2 V > 15 GeV. In the H → 4 decay, we model the C HW B = 1 contribution to the SM, as shown in Fig. 2 with cross section decomposition presented in Table IV. In VBF or V H, we model the C H W B = 10 or C HB = 100 contribution to the SM, as shown in Fig. 3 or Fig. 4, with the cross section decomposition presented in Table V or Table VI. The size of anomalous contributions is chosen to be large compared to SM for visibility of their contributions.
In the H → 4 process, the larger and the smaller invariant masses of the dilepton pairs m 1 and m 2 are the two observables representing q 2 1 and q 2 2 . In Fig. 2, there are clear peaks towards m 2 → 0 in the case of couplings with photons, HZγ and Hγγ, a feature to which we will return in Sections IV, V and VI. In the case of Hγγ, this extends to m 1 → 0 as well. Modeling such contributions becomes essential, and we will discuss extensions of such modeling to m < 1 GeV later. Moreover, in analysis of experimental data, detector effects change significantly for either γ * or Z intermediate states, and dedicated simulation of such effects with the full detector modeling becomes important. In Fig. 2, the m 1 and m 2 distributions are shown separately in the H → 4e/4µ and H → 2e2µ decays. The interference of two diagrams with permutation of identical leptons in the case of H → 4e/4µ leads to suppression of the peaks at m 1 → 0 and m 2 → 0 in the case of Hγγ. This feature becomes important in analysis of the Hγγ couplings.
Since Zf f couplings have been constrained with precision EW data, we do not allow their change in these studies and assume that modification of other operators, not contributing to the H boson couplings, can compensate any possible shift of the Zf f couplings due to C HW B . However, in Fig. 2 we also show distributions with modification of the Zf f couplings, indicated with δg Zf f . Corrections to the multidimensional angular distributions are expected due to non-zero values of the R i and A f parameters discussed in Refs. [19,23,28]. These corrections are visible in the projection on the Φ observable in Fig. 2, but are very small for any practical purpose with the typical values of C HW B in the present studies. These corrections become sizable with larger values of C HW B .
In the VBF process, we can calculate the q VBF 1,2 = −q 2 1,2 values using the momenta of the fully reconstructed H boson and two jets and using the direction of incoming partons along the proton beams. In Fig. 3, there is a clear preference of lower q VBF 1,2 values in the case of couplings with photons. There is a strong correlation between the q VBF 1,2 values and the transverse momentum p T of the jets, which leads to different detector effects. We note the asymmetric distribution of the Φ VBF angle in Fig. 3, which is most visible in the HZZ process but can also be seen in the combined distribution. This happens due to interference of the CP -even SM amplitude and CP -odd C H W B = 10 contributions.
In the V H process, q 2 1 and q 2 2 represent the V H and the V → + − invariant masses, respectively. There are particularly dramatic effects in the m distribution, shown in Fig. 4, where the virtual photon γ * results in the  low-mass enhancement, as opposed to the peak at m Z . A dedicated analysis of the small invariant masses in the γ * H production may be needed for effective EFT analysis of the process. Another approach to study anomalous H boson couplings involving photons is analysis of the γH process, which distinguishing feature is a high-momentum on-shell photon associated with the H boson. In the LO topology, where the γH system has no transverse boost, the transverse momentum of either photon γ or the H boson is a dependent observable and the three primary measurements are the rapidity y and the invariant mass m γH of the m γH system, and the angle θ 1 formed by the outgoing photon with respect to the direction of incoming quark in the γH rest frame. This angle is also defined in Fig. 1, where V 2 = γ and which does not have a subsequent decay. While it is not possible to distinguish the incoming quark and antiquark on event-by-event basis, on average the boost direction of the γH provides the preferred direction of the quark, and we use this to define θ 1 . However, determination of the cos θ 1 sign becomes important only in the special case of the forward-backward asymmetry discussed below. The ability to determine the cos θ 1 sign is a function of y and has been discussed earlier [45,91].
In Fig. 5, the m γH and cos θ V H 1 distributions are shown for the g Zγ 2 , g γγ 2 , g Zγ 4 , or g γγ 4 anomalous couplings, and for the mixture of g Zγ 2 and g Zγ 4 contributions with a complex phase of g Zγ 4 . The m γH distributions differ somewhat between the HZγ and Hγγ couplings, due to the difference between the intermediate Z * and γ * . The cos θ V H 1 distributions follow the 1 + cos 2 θ 1 expectation for all real couplings. This expectation can be traced back to Eq. (A2) in Ref. [28], where A 00 = 0, which must be averaged over Φ and cos θ 2 . The situation becomes similar to the angular distribution in the H → Zγ decay, described with the same angular parameterization, as shown in Fig. 15 of Ref. [28], and where the forward-backward asymmetry may be generated with the mixture of CP -odd and CP -even couplings and in the presence of a complex phase.
The size of the forward-backward asymmetry is proportional to the A f parameter defined in Ref. [28] for Zf f  couplings, which is 0.15 for the lepton couplings, but is as large as 0.67 and 0.94 for the up and down type quarks. Therefore, despite the sizable dilution in the measurement of the cos θ 1 sign, the forward-backward asymmetry is strongly pronounced in Fig. 5 and could be measured in experiment once γH production is observed. The A f parameter is zero for the photon couplings γf f , and such an effect is not possible in the mixture of couplings involving g γγ 2 and g γγ 4 . Since non-trivial forward-backward asymmetry appears only in the special case of complex couplings, we do not consider this asymmetry further in this work, but we point out that such a study is in principle possible.

IV. PARAMETERIZATION OF CROSS SECTIONS
In this Section, we discuss the relationship between the coupling constants and the cross section of a process involving the H boson. In Ref. [45], we calculated the scaling factors for the partial decay widths in the nine dominant H boson decay modes as a function of anomalous couplings a i , including the H → gg, γγ, and Zγ decays, by resolving the loop contributions. However, we omitted point-like contributions such as g γγ 2,4 and g Zγ 2,4 due to their relatively lower importance in the VBF, V H, and H → 4f processes. Such couplings could be generated by a heavy quark Q with mass m Q M H . We assume that its couplings to the H boson are κ Q andκ Q , the number of colors N c , the electric charge Q, and the weak isospin projection T 3L . This special model allows us to derive the point-like interactions and relate those to the partial decay widths. While derivation applies to this special case, the final expression in terms of the g gg 2,4 , g γγ 2,4 , and g Zγ 2,4 couplings becomes generic and remains valid for any new physics in the loop, generated by any combination of heavy fermions or bosons. Therefore, the resulting expressions are applicable to the general treatment of these loops in the EFT approach.
First, we recall that in the narrow-width approximation for on-shell H boson production and decay, the cross section can be expressed as where the total width Γ tot = Γ known + Γ other representing decays to known particles and other unknown final states, either invisible or undetected in experiment. In the following we will focus on decay to the known SM particles which can be expressed as a sum of all partial decay widths as where R f is the scaling factor as function of the coupling constants a i , and Γ SM f is the SM value of the partial decay width in the final state f .  In the following, we rely on JHUGen framework implementation, discussed in Section III and Ref. [45], to derive the loop contributions of the SM particles and the heavy quark Q to the scaling factor R gg , for both CP -even and CP -odd couplings. The CP -even coupling contributions of the quarks and W boson to R γγ and R Zγ are derived with HDECAY [92]. The CP -odd contributions to R γγ are calculated with the JHUGen framework in a manner analogous to R gg . The CP -odd contributions to R Zγ are calculated using cHDECAY [93].
The ratio of the decay width to the SM expectation in the H → gg process [45] is found to be The κ Q andκ Q couplings are connected to the g gg 2 and g gg 4 point-like interactions introduced in Eq. (1) through One can rewrite Eq. (13) in terms of the g gg 2 and g gg 4 couplings in place of N c κ Q and N cκQ by substituting Eq. (14). Even though Eq. (14) is derived in the special case of a heavy quark, the resulting expression of R gg as a function of g gg 2 and g gg 4 and other terms is valid for any heavy particles in the loop that generate these point-like interactions. The latter observation allows us to obtain the value of the effective g gg 2 coupling which leads to the SM cross section in the gluons fusion process. By setting all couplings, other than g gg 2 , to zero and R gg = 1 in Eq. (13), we obtain The g gg,SM 2 value differs by only 1.5% from the value that one would obtain in the heavy top mass limit by setting κ Q = 1 and N c = 3 in Eq. (14), and the sign follows the prediction in this limit.   15), substitute κ Q andκ Q for g gg 2 and g gg 4 , and obtain For the H → γγ final states, we include the W boson in addition to the top, bottom, and heavy Q quarks in the loop and obtain 2 For the contribution of a heavy quark in the loop we find Following the idea described above for R gg , one can rewrite Eq. (17) in terms of the g γγ 2 and g γγ 4 couplings in place of N c Q 2 κ Q and N c Q 2κ Q by substituting Eq. (18). The final expression of R γγ as a function of g γγ 2 and g γγ 4 and other terms is again valid for any heavy particles in the loop, fermions or bosons, that generate these point-like interactions. By setting all couplings other than g γγ 2 to zero and R γγ = 1 in Eq. (17), we obtain the effective coupling which leads to the SM cross section The g γγ,SM 2 value differs slightly from 0.00400 obtained from the general expression of the SM loops derived from Refs. [94,95] and shown in Eq. (20). The difference could be explained by the higher-order effects incorporated in Eq. (17) and the fact that in our approach we match the SM rate R γγ = 1. The sign in Eq. (19) follows Eq. (20).
where the one-loop functions are given by and An approximate way to express Eq. (17) with point-like interactions only would be to follow the idea used to create Eq. (16) and substitute the SM couplings with g γγ,SM 2 from Eq. (19), substitute κ Q andκ Q for g γγ 2 and g γγ 4 , and obtain For the H → Zγ final states, for the coupling of the heavy Q quark to the Z boson, we introduce the following parameter which corresponds to the following values for the SM parameters of the top (T 3L t = +1/2, Q t = +2/3) and bottom We obtain For the contribution of a heavy forth generation quarks in the loop we find We note that the effective value of g Zγ 2 for a heavy quark Q which reproduces the SM partial width, is The g Zγ,SM 2 value differs slightly from 0.00724 obtained from the general expression of the SM loops derived from Refs. [94,95] 3 and shown in Eq. (29). As before, the difference could be explained by the higher-order effects incorporated in Eq. (26) and the fact that in our approach we match the SM rate R Zγ = 1. The sign in Eq. (28) follows Eq. (29).
where the one-loop functions are given by A Zγ 35. An approximate way to express Eq. (26) with point-like interactions only would be to substitute the SM contributions with an effective coupling g Zγ,SM 2 from Eq. (28), substitute κ Q andκ Q for g Zγ 2 and g Zγ 4 , and obtain In the above calculation, the H → γ * γ process is not included, for which the full loop calculation with anomalous couplings is not available. For the H → ZZ/Zγ * /γ * γ * → four-fermion final state, the full one-loop calculation with anomalous couplings is not available either. The EW loop corrections under the SM assumption are discussed in Section V. A more careful treatment of the singularities appearing in the presence of anomalous couplings in both of the above cases is discussed in Section VI. For the leading tree-level contributions, we derived the R ZZ/Zγ * /γ * γ * parameterization in Ref. [45], in which case we set g Zγ 2 = g Zγ 4 = g γγ 2 = g γγ 4 = 0 to avoid collinear singularities. In the following, we introduce these four couplings and avoid singularities in the γ * → 2f transition with the finite fermion mass threshold (1) and rely on the κ Zγ 2 and κ ZZ 1 = κ ZZ 2 parameters to express the scaling factor as 4 Equation (31) covers all final states with Z/γ * → qq and + − with quarks and charged leptons, while neutrinos are included with Z → νν. The treatment of qq hadronization with the low-mass resonances is not included here, and is discussed in more detail in Section VI. The interference between the CP -odd and CP -even contribution integrates out to zero, as reflected in the zero terms in Eq. (31).
Let us conclude this Section by discussing the cross section of the qq → γH process as a function of the anomalous couplings summarized in Table VIII. Detecting or setting limits on this process will be of interest for constraining the 4 There is a sign change of the κ Zγ 2 coupling when compared to coefficients in common with Ref. [45], because here we use the convention Bµ, as discussed in Section II. following couplings, as discussed in Section III: where the reference cross section is σ γH ref = 1.33 × 10 4 fb. We will investigate this channel further in Section VII.

V. LOOP-INDUCED STANDARD MODEL CONTRIBUTIONS
The NLO EW corrections from the Prophecy4f and HAWK generators have been widely used in calculations of the H boson production and decay cross sections at the LHC and included in the LHC Higgs Working Group recommendations [39]. The corrections are generally positive in the H → 4 process [96] and negative in the VBF and V H processes [97]. Differential distributions also show growth of these effects at higher energy, such as at high transverse momentum p H T of the H boson in the case of VBF and V H production, as expected for the well-known EW Sudakov enhancement. Our goal here is to reexamine some of these effects, focus on certain kinematic distributions, and compare the NLO EW effects to those generated by the EFT operators. In particular, we also produce kinematic distributions with JHUGen at LO, and introduce effective g γγ,SM 2 and g Zγ,SM 2 couplings to model what one can call pseudo-EW corrections. Both Prophecy4f and HAWK include the interference of the loop-induced contributions with the Born process as dictated at NLO accuracy, but do not include squared contributions, which are formally of higher order. Nonetheless, these squared terms may be comparable to or larger than the interference contributions, and we examine this with the effective g γγ,SM be found in Eqs. (15), (19), and (28), respectively. As we illustrate below, these couplings are inadequate for modeling loop effects in decays to virtual vector bosons, such as H → γ * γ/Zγ → ff γ or H → γ * γ * /Zγ * /ZZ → ff f f . The non-trivial q 2 dependence of the effective HV V vertex cannot be described this way and the g ZZ 2 and other tensor structures appearing in Eq. (1) are not represented. Similar considerations apply to the VBF and V H production, and also to the γH production, as discussed in Section III. To make these statements in a quantitative way, we compare a simluation of these effective couplings to a complete modeling of the NLO EW effects.
In order to model the SM loop corrections in the H → ff f f process, we employ the Prophecy4f generator, and in VBF and ZH production we use HAWK. In both cases, the NLO EW corrections can be applied to the process of interest and compared to the LO simulation. We note that HAWK provides all results in the form of binned distributions, since unweighted events are not available and events are not stored in LHE format [98]. This may complicate analysis and comparison of generated events, since different code would have to be employed in calculating the observables. Therefore, we have introduced a software interface which writes weighted events from HAWK simulation in the LHE format, and all further analysis is performed in a unified way. Moreover, photon bremsstrahlung leads to smearing of kinematic distributions. In order to disentangle photon radiation from purely EW effects in kinematic distributions in both the H → 4 and V H with V → + − processes, we have introduced a recombination algorithm in analysis of events written in LHE format. The four-momentum of the associated photon is added to the nearest lepton in this algorithm. We note that NLO QCD+EW predictions for ZH production have been recently implemented in the POWHEG framework [99], but are not included in this study.
We start with the study of NLO EW corrections in the H → 2e2µ decay process. In Fig. 6, the LO and NLO EW modeling of the process is shown, as generated with the Prophecy4f generator, and compared to ad-hoc loop correction with g γγ,SM 2 and g Zγ,SM 2 with JHUGen. The overall correction to the decay width is +1.5%, as shown in Table II. The size of the effect with the pseudo-EW correction of −0.6% is similar, but does not reproduce the sign. However, including the quadratic terms with the pseudo-EW corrections appears important, as there is a growing importance of these effects at q 2 → 0, which increases the correction by +2.6%. The effect of linear terms appearing with the proper NLO EW corrections is most pronounced in the intermediate m 1 and m 2 ranges, away from the pole of on-shell Z. This is where the effect of interference is most pronounced. Overall, we conclude that the pseudo-EW corrections only roughly model the effect in the H → 2e2µ process of an order of magnitude, but are not adequate to describe the proper EW corrections. Nonetheless, they also indicate that the quadratic terms may become sizable and more important than the linear terms at the values of m 2 of a few GeV or below.
The study of NLO EW corrections in the VBF process is shown in Fig. 7 and in the V H process in Fig. 8, where the LO and NLO EW modeling of the process is shown with the HAWK generator, and compared to ad-hoc loop correction with g γγ,SM 2 and g Zγ,SM 2 with JHUGen. The selection requirements are similar to those in Section III, except that in VBF we do not place a requirement on q 2 V and in V H we apply a looser requirement m > 0.1 GeV and a tighter requirement p T > 5 GeV. The overall NLO EW correction is negative in the range of 6 − 7%, as shown in Table II  about −1.2% correction in the V H process and small growth of the effect with p H T , but no sizable effect in the VBF process. In both VBF and V H, there is no evidence of importance of the quadratic terms with the g SM 2 expansion. We conclude that the pseudo-EW corrections are not adequate in the VBF and V H processes, even if in the V H process cross section modifications appear in the same direction.
Given the inadequacy of the pseudo-EW corrections, we have investigated an approximate approach of re-weighting the LO EW simulation with a dynamic k factor, which is a ratio of the NLO and LO EW kinematic distributions. This approach can capture the main features of the correction, such as its growth with p H T , if this quantity is used as the kinematic distribution in re-weighting. However, this approach does not guarantee adequate modeling of the other kinematic distributions simultaneously, if those are not also used in re-weighting. For example, we found that re-weighting based on p H T does not bring angular distributions to an agreement. The importance of the NLO EW corrections will become evident in the actual analysis of LHC data. The precision of existing LHC constraints [3,[46][47][48][49][50][51][52][53][54][55][56][57][58][59][60] is not sufficient for reaching the NLO EW effects appearing in the SM. Therefore, accurate modeling of such effects may not appear as critical at present. However, with growing precision of experimental measurements a careful investigation of NLO EW effects on kinematic distributions will become important.

VI. CALCULATION OF THE WIDTH IN THE PRESENCE OF INTERMEDIATE PHOTONS
The presence of intermediate photons in H boson decays to fermions via γ * → f + f − can lead to sharp peaks in the spectrum when q 2 = (p f + + p f − ) 2 is small. For decays into leptons, these peaks are cut off by the physical lepton masses. In the case of quarks, non-perturbative effects wash out this peak structure and introduce hadronic resonances instead. In the following, we introduce a procedure, based on matching amplitudes in the collinear limit, to handle these singularities in a way which allows their efficient numerical evaluation. We draw parts of this description from Ref. [62]. Technical details can be found in Appendix B.
We write the partial H boson decay widths as where we introduced the symbols (j pp ) = j δ pp for symmetry factors for identical particles. Explicit parameterization of the phase spaces, using appropriate variables, yields and where simply dPS (2) (q 2 , m 2 1 , m 2 2 ) = dcos θ dφ/(2q 2 ). We split the invariant mass integrations in Eqs. (36,37) into a low and high virtuality region separated by an arbitary parameter µ 2 M 2 . In this form, we can apply the collinar approximation to the squared in the low virtuality region and analytically integrate over q 2 γ [62]. The high virtuality region does not contain sharp peaks and can be treated numerically in a standard manner. As a result, the partial decay width for H → 2 γ can be written as where the left-hand side is independent of µ 2 . We note that the low virtuality region is conveniently expressed in terms of an analytic function containing a potentially large log(m 2 ) times the width Γ H→γγ , which can be straight-forwardly obtained using numerical methods for any combination of anomalous couplings (cfg. Eq. (17)). Contributions from low virtuality Z bosons decaying into 2 are parametrically suppressed by µ 2 M 2 H . In a similar fashion, we obtain the results for the H boson decays into four leptons If the H boson decay occurs into quark final states, the low virtuality region is affected by sizable non-perturbative effects, which can be related to the experimentally measured quantity ∆α (5) had (M 2 Z ) via a dispersion relation and unitarity [62]. The value of ∆α (5) had (M 2 Z ) = (276.11 ± 1.11) × 10 −4 [100,101] has been extracted from the low energy region of the ratio σ(e + e − → hadrons) σ(e + e − → µ + µ − ) and is related to the fine structure constant via α(s) = α(0) (1 − ∆α(s)). Summing over the five light quark flavors, labelling q Γ H→2qγ = Γ H→2jγ , and choosing and similar for Γ H→4j . In order to illustrate the performance of Eq. (39), we model the H → 2 γ decay with the HV V couplings corresponding to g γγ,SM 2 and g Zγ,SM 2 as defined in Eqs. (19) and (28). We scan the value of cutoff µ 2 in Eq. (39) from the threshold value of (2m ) 2 up to (5 GeV) 2 for both electrons ( = e) and muons ( = µ). Figure 9 shows the first two terms appearing in Eq. (39), the value of Γ H→2 γ with the requirement q 2 2 ≥ µ 2 and the value of Γ H→γγ multiplied by the µ 2 -dependent factor. The other terms in Eq. (39) can be neglected in this comparison. With the µ 2 increasing, the former cross section falls while the latter rises, leading to a constant cross section of the H → 2 γ process, as one would expect to observe for the proper modeling of the effect.  41) with the couplings discussed in text. The first term in each equation is shown with the requirement q 2 2 ≥ µ 2 common for electrons and muons (blue diamonds, ). The following terms multiplied by the µ 2 -dependent factors are also shown: the Γ H→γγ (red circles, •), the Γ H→2eγ with the requirement q 2 2e ≥ µ 2 (magenta triangles pointing down, ), and the Γ H→2µγ with the requirement q 2 2µ ≥ µ 2 (brown triangles pointing up, ). The sum of the cross sections is also shown (black squares, ).
In order to illustrate the performance of Eqs. (40) and (41) , respectively, as used in the previous test. Both in this test and in the previous test of Eq. (39), the rationale is to model the processes involving virtual photons which exhibit growth with q 2 → 0. We choose to illustrate the performance with CP -even couplings in one case and CP -odd couplings in the other, but the procedure has been validated to work in any combination of couplings. The κ Zγ 2 coupling formally involves a virtual photon. However, due to appearance of q 2 γ in the numerator of the tensor structure in Eq. (1), this coupling does not exhibit divergence with q 2 → 0. Its behavior is similar to the terms with the virtual Z and can be absorbed in the corresponding terms.
As before, we scan the value of cutoff µ 2 in Eq. (40) for the H → 2e2µ process and in Eq. (41) for the H → 4e and 4µ processes from the threshold value of (2m ) 2 up to (5 GeV) 2 . The threshold value is defined for muons in the case of H → 4µ and electrons in the other two cases. Figure 10 shows the partial decay width of the process H → 2e2µ calculated using Eq. (40), and H → 4e and 4µ calculated using Eq. (41). Several terms in the corresponding equations are isolated: the first term in each equation is shown with the requirement q 2 2 ≥ µ 2 common for electrons and muons; the Γ H→γγ multiplied by the µ 2 -dependent factor, the Γ H→2eγ with the requirement q 2 2e ≥ µ 2 , and the Γ H→2µγ with the requirement q 2 2µ ≥ µ 2 , where both Γ H→2 γ are multiplied by the µ 2 -dependent factor as well. With the µ 2 increasing, the first term falls while the other terms generally rise. Again, we obtain a constant cross section of the four-lepton process, as one would expect to observe for the proper modeling of the effect.
The combination of formulas in Eqs. (39) and 42 for Γ H→2f γ and Eqs. (40), (41), and similar ones involving hadronic jets for Γ H→4f should allow one to handle low-q 2 singularities and hadronic structure with efficient numerical evaluation.

VII. EXPECTED CONSTRAINTS ON THE COUPLINGS WITH PHOTONS
We continue by investigating the on-shell production and decay of the H boson with its couplings to weak vector bosons in the VBF, V H, and H → V V → 4 processes. There has already been extensive study of the HV V couplings, both by experimental collaborations and in phenomenological work. However, there was no conclusive study on the effects of the photon contribution in the production topology. In our previous work in Ref. [45], we pointed out that the decays H → γγ and Zγ with on-shell photons provide constraints on the g γγ 2 , g γγ 4 , g Zγ 2 , and g Zγ 4 couplings which are stronger than those that can be obtained from the VBF, V H, and H → 4 processes. For this reason, the analysis of multiple operators was simplified by setting those four couplings to zero. The constraints from the decays with on-shell photons can be illustrated with the simplified partial decay width expressions in Eqs. (23,30). This effect in the H → 4 was studied with LHC data [48] and with phenomenological tools [25,35]. In the following we re-examine the H → 4 decay and investigate the VBF and V H processes.
First, we would like to point to the effect already observed in Tables IV, V, and VI. Let us focus on any of the three operators C HW , C HW B , or C HB where the g ZZ 2 , g Zγ 2 , and g γγ 2 couplings contribute, or equivalently any of C H W , C H W B , or C H B , where g ZZ 4 , g Zγ 4 , and g γγ 4 contribute. While the g ZZ 2,4 contributions to the VBF and V H processes are comparable and sometimes even dominant for a given C HX operator, their contributions to the H → 4 decay appear to be negligible in comparison. For photon couplings, the reverse is the case: their contribution to decay is much larger than to production. Therefore, the photon couplings have relatively higher importance in decay compared to the production processes. This still does not tell us if the photon couplings are better constrained in production or decay, and this is what we investigate below.
In the following, prospects with either 3000 fb −1 (HL-LHC) or 300 fb −1 (full LHC) are studied at a 13 TeV collision energy. We use the JHUGen simulation to model the VBF, W H, ZH/γ * H, γH, and gluon fusion production with the decay H → 4 . We include the effective background with the POWHEG [102] simulation of the qq → 4 process, which we scale to match the contributions of other processes as found in experiment, such as gg → 4 and Drell-Yan Z production [103]. The detector effects are modeled with ad-hoc acceptance selection and empirical efficiency corrections, and the lepton and hadronic jet momenta are smeared to achieve realistic resolution effects. The following selection requirements are applied: We target the optimal analysis of four anomalous couplings expressed through the fractional contributions to the H → 2e2µ process f Zγ g2 , f γγ g2 , f Zγ g4 , and f γγ g4 , with the approach similar to Ref. [45]. The cross section fractions are defined following where the α (f ) nn coefficients are introduced in Eq. (11) for the final state (f ). The numerical values of these coefficients are given in Table III, where they are normalized with respect to the α 11 coefficient, corresponding to the cross section calculated for g 1 = 1. The α nn are the cross sections for g n = 1. In Table III, we also quote the cross section ratios defined in VBF and V H production for comparison. As noted above and evident from this table, the relative importance of photon couplings is higher in decay compared to production.
For simplicity, we set the contributions of the other anomalous contributions to zero: f ZZ g2 = f ZZ g4 = f ZZ Λ1 = f Zγ Λ1 = 0. This assumption will provide tighter constraints than one could achieve otherwise, but this is sufficient for comparison to the precision obtained from H → γγ and Zγ following Eqs. (23) and (30). We use the relationship with ∆M W = 0 in Eq. (5) for the SM-like contribution, but we do not enforce the SU (2) × U (1) symmetry in Eqs. (6-9) and instead set g W W This is done to isolate the contributions with genuine photon couplings, which exhibit the features of virtual photons. These are the contributions which also appear and are constrained in the H → Zγ and γγ processes.

A. Expected constraints on photon couplings from H → V V → 4 decay
We build the analysis following the MELA approach with the two types of optimal discriminants, using the full kinematic information in the four-lepton decay. There are four discriminants D Zγ g2 , D γγ g2 , D Zγ g4 , and D γγ g4 which are designed to separate the pure anomalous contributions from the SM-like, and there are four interference discriminants D Zγ int , D γγ int , D Zγ CP , D γγ CP which isolate interference of the SM with the same four anomalous contributions. The full available information is used in calculating the discriminants, and further details on the MELA approach can be found in Ref. [45] and references therein. In this and further analyses discussed below, events are split into the H → 2e2µ, , and f Zγ g4 using associated production and H → 4 decay with 3000 (300) fb −1 data at 13 TeV. Two scenarios are shown: using MELA observables with production and decay (red) or production only (blue) information. The dashed horizontal lines show the 68 and 95% confidence level (CL) regions.
4e, and 4µ categories, which is an important aspect because the relative fractions of events between these types change with anomalous couplings, due to interference of diagrams with identical leptons in the 4e and 4µ final states. In the end, for each event in a category j, a set of observables x is defined. Since analysis of decay information is essentially independent from the production mechanism, we model kinematic distributions using simulation of the gluon fusion process. However, in Section VII B we will also show how the full production and decay information can be taken into account. The probability density function for the signal decay process in gluon fusion production, before proper normalization, is defined as ggH : where x are the observables and f gn are the cross-section fractions of the couplings, K = 5 for the four anomalous couplings and one SM coupling. Equation (44) is obtained from Eq. (11), where the width and Hgg couplings are absorbed into the overall normalization. In gluon fusion production, the electroweak HV V couplings appear only in decay. Therefore, there are 15 terms in Eq. (44). To populate the probability distributions, we use a simulation of unweighted events of several samples that adequately cover the phase space and re-weight those samples using the MELA package to parameterize the other terms. As in Ref. [45], we implement a cutting planes algorithm [104] using the Hom4PS [105][106][107] and Gurobi [108] programs to ensure that the probability density function remains positive definite for all possible values of f gn . With 3000 fb −1 data at 13 TeV, we expect about 4500 events reconstructed in the H → 4 channel. In Fig. 11 we show the expected constraints on the four parameters of interest expressed through effective factions f γγ g2 , f γγ g4 , f Zγ g2 , and f Zγ g4 . However, before discussing the results, we should turn to analysis of production information.

B. Expected constraints on photon couplings from VBF and V H production
While analysis of the H → 4 decay can be performed inclusively, the VBF and V H production channels require analysis of associated particles. Therefore, we follow the approach adopted in our earlier studies [45] to categorize events using the jet information to create the two-jet categories of events enhanced in either VBF or V H events with hadronic decay of the V , respectively. In addition, the leptonic V H, one-jet VBF, and boosted categories are defined, where p 4 T > 120 GeV is used to select the latter [45]. The remaining events constitute the untagged category. We build the analysis following the MELA approach with the two types of optimal discriminants, as discussed for decay above, using the full kinematic information in both the production and the four-lepton decay. In the VBF and V H topologies with two associated jets, both production and decay information are used, except for the interference discriminants, where production information is chosen. In the untagged category, the H → 4 information is used as in Section VII A, and in the other categories the transverse momentum p 4 T is used. In the case of the VBF and V H processes, the HV V coupling appears on both the production and decay sides. Therefore, the amplitude squared has a product of four couplings. Equation (45)  using decay and production information as shown in Fig. 11. 70 terms: where the notation follows Eq. (44). In addition to the joint analysis of production and decay, we also perform an analysis with production information only. In order to achieve this, no decay information is used in the construction of discriminants, and only the total yield of events is used in the untagged category. It is not possible to completely decouple the analysis from decay information, as for example the relative fraction of 2e2µ events is sensitive to couplings in the decay amplitude. For example, there is a strong destructive interference between the two diagrams with permutation of identical leptons in the H → 4e/4µ decay with the Hγγ couplings, as illustrated in Fig. 2 and Section III, which leads to a modification of the 2e2µ yield faction. However, such dependence on anomalous couplings is greatly reduced with consideration of yields only.
In Fig. 11 we show the expected constraints on the four parameters of interest f γγ g2 , f γγ g4 , f Zγ g2 , and f Zγ g4 with 300 and 3000 fb −1 data at 13 TeV, where for comparison constraints from production information alone are shown separately from the full constraints using both decay and production. The information contained in decay can be deduced from the difference of the two constraints, with one exception. The analysis of production information is sensitive to the relative fraction of the H → 2e2µ, 4e, and 4µ events. We observe that in the case of the H → γ * γ * → 4e and 4µ processes, there is a sizable effect of interference between diagrams with permutations of the leptons. This effect in decay competes with information from production in constraining g γγ 2 and g γγ 4 . We do not find such an effect to be important in constraining g Zγ 2 and g Zγ 4 . What we observe is that the constraints from decay are significantly more powerful than from production when both are analyzed with the same channel H → 4 . This happens for two reasons. First of all, all reconstructed events contain decay information, while only a small fraction carry production information. Second, the ratio α (i) nn /α 11 is reduced compared to the same ratio in decay α (f ) nn /α 11 for the photon couplings. This effect is in contrast to the trend observed for the ZZ couplings (see, for example, Fig. 10 of Ref. [45]), as indicated with the g ZZ 4 example in Table III. This trend explains the tighter constraints on the anomalous ZZ couplings using production information, as opposed to decay. For the same reason, the constraints on the photon couplings (Hγγ and HZγ) are tighter in decay compared to production.
In Fig. 12, these results for the joint analysis of production and decay are interpreted in terms of constraints on the g i couplings. As indicated above, these constraints are dominated by decay information, and the results would look similar if only decay information were employed. This interpretation uses the full expression in Eq. (11), including the total H boson width dependence on anomalous couplings appearing in the denominator, using expressions obtained in Section IV. It is assumed that there no decays of the H boson to unknown particles.
We should point out that while decay information is limited to H → V V channels only, production information can be obtained by combining all possible decay channels of the H boson, such as H → 4 , γγ, bb, τ + τ − , and W + W − . In this study, we investigate only the H → 4 channel, but the relative importance of production information will  using the γH production channel.
increase as other channels are analyzed. This observation is also valid for analysis of the γH production, discussed next. At the same time, the H → 4 channel can also be further optimized for the measurements of the photon couplings by relaxing the invariant mass and transverse momentum constraints on the leptons. For example, the requirement |m | > 12 GeV can be relaxed to |m | > 4 GeV, or even further, with significant gain in sensitivity to the Hγγ and HZγ couplings, as can be seen from the invariant mass m 1 and m 2 distributions in Fig. 2.
C. Expected constraints on photon couplings from γH production Setting constraints on or measuring the rate of the γH production is interesting on its own, as this production mechanism of the H boson has not been tested on the LHC, in part because its SM rate is not accessible yet. In addition, we would like to assess the feasibility of the g γγ 2 , g γγ 4 , g Zγ 2 , and g Zγ 4 coupling measurement in this production process. In order to simplify these estimates and because NLO EW event simulation is not available for the γH process yet, we assume that the SM contributions are small compared to the accessible values and thus can be absorbed into these effective point-like couplings. The validity of these assumptions can be checked against the expected constraints. Single loop calculations of the SM production cross section of qq → γH predict a cross section of less than 5 fb at the LHC [109], which corresponds to less than two H → 4 event at the HL-LHC before any selection requirements are applied.
Experimentally, the main distinguishing feature of the γH production mechanism is a high-momentum isolated photon, and we found that the requirement p γ T > 400 GeV keeps about 12% of signal events, which corresponds to about 0.9 signal events for g γγ 2 = 0.1 and 0.11 background events with H → 4 and 3000 fb −1 data at 13 TeV. Therefore, for feasibility studies, we identify an additional category of events with a good isolated photon candidate associated with the H boson candidate and the above p γ T requirement. Since kinematic features of events do not differ between CP -odd and CP -even couplings and the difference between the HZγ and Hγγ couplings is weak, we perform a simple fit for excess of events over background without using kinematic distributions. The expected number of signal events is expressed through couplings as the product of coupling modifies on the decay and production sides. These expressions are similar to those in Eqs. (31) and (32), but differ in two aspects. Equation (31) has to be adjusted for the H → 4 final state instead of the inclusive four-fermion final state and include the effect of acceptance requirements on the lepton quantities, as listed above. Equation (32) has to take into account the effects of the photon acceptance efficiency.
With the above assumptions, the expected two-dimensional constraints on (g γγ 2 , g γγ 4 ), (g Zγ 2 , g Zγ 4 ), and (g Zγ 2 , g γγ 2 ) are shown in Fig. 13. The expected γH constraints with the H → 4 channel alone are not as powerful as those obtained from decay and shown in Fig. 12. However, while this difference is sizable in the case of the Hγγ couplings, the difference is not as large in the case of the HZγ couplings. Moreover, these constraints are comparable and even better than those obtained from production information in the VBF and V H channels. In both cases, significant gain will result from the analysis of the other H boson decay channels, which are not considered in this feasibility study. Therefore, it is important to proceed with analysis of VBF, V H, and γH production in all accessible H boson final states.
Searches for heavy resonances decaying into a photon and a hadronically decaying H boson have been performed on   (23) and (30) and the HL-LHC projection for experimental measurements of R γγ and R Zγ from Ref. [112]. Inclusion of the H → 4 data with decay and production information helps in resolving degenerate solutions on a ring of a given radius.
the LHC [110,111], where the topology of the final state overlaps with the process of our interest, but the interpretation of the results is very different. A reinterpretation of these resonance-search results in terms of the EFT couplings of the H boson was attempted in Ref. [84]. However, this study makes crude approximations in its interpretation of LHC data, uses a different set of parameters in a different basis, applies additional symmetry considerations, and is limited to CP -even couplings only, making a direct comparison with our results difficult.

D. Expected constraints on photon couplings from decays with on-shell photons
The above results can be compared to possible constraints from the measurements of R γγ and R Zγ , defined in Section IV. The projection of experimental measurements of the H boson branching fractions to 3000 fb −1 of LHC data has been performed in Ref. [112]. In particular, Table 37 of this reference estimates R γγ 1.00 ± 0.05 and R Zγ 1.00 ± 0.24 at 68% CL. However, these measurements are estimated using the coupling modifier framework (κframework), where the tensor structures of interactions of the H boson and the kinematic distributions are assumed to be the same as in the SM. Though the kinematic distributions in decays H → γγ and Zγ are not affected by anomalous couplings, other aspects of the analyses, such as distributions of associated particles, would be affected. Nonetheless, one can use these estimates as optimistic expectations of constraints on R γγ and R Zγ with anomalous contributions.
Relating the experimental constraints on R γγ and R Zγ to constraints on g γγ 2 , g γγ 4 , g Zγ 2 , and g Zγ 4 is not trivial, as various loop contributions are possible, as indicated in Eqs. (17) and (26). In particular, modifications of the H boson couplings to fermions and W boson, which are the dominant SM contributions to the loop, cannot be disentangled from the effective point-like couplings indicated with the κ Q andκ Q couplings. The latter are related to g γγ,Zγ 2,4 with Q and R Q in Eqs. (18) and (27). However, assuming that the H boson couplings to fermions and W boson can be constrained to high precision from other measurements, one can set those to SM values for the purpose of this comparison. In such a case, R γγ and R Zγ could be expressed directly through the g γγ 2,4 and g Zγ 2,4 couplings, respectively, without further complication, as shown in Eqs. (23) and (30).
With the above assumptions, the expected two-dimensional constraints on (g γγ 2 , g γγ 4 ) and on (g Zγ 2 , g Zγ 4 ) are circles in the 2D plane of the two couplings, as follows from Eqs. (23) and (30), respectively. All points on the circle of a given radius are equally likely, because one cannot distinguish between the CP -odd and CP -even couplings from the rate information alone. However, the addition of H → 4 decay and production information, discussed above, will help to resolve degenerate solutions, as shown in Fig. 14. The improvement from the inclusion of the H → 4 data is more visible in the case of the Hγγ couplings. As we would expect, constraints from the H → γγ and Zγ rates are more restrictive than the constraints obtained from either decay or production information shown in Figs. 12 and 13.
As we reach precision on anomalous Hγγ and HZγ couplings from the H → 4 decay and from VBF, V H, and γH production comparable to that from the H → γγ and Zγ decays, resolving the loop contributions in these production and decay processes will become important. This is similar to the discussion of H → γγ and Zγ in Section IV, but taking into account non-trivial q 2 -dependence affecting kinematic distributions. This is equivalent to the NLO EW corrections discussed in Section V in application to the SM processes, but would require consideration of anomalous couplings.

VIII. SUMMARY
We have presented a study of electroweak production of the H boson in VBF and V H and its decay to two vector bosons, with a focus on the treatment of virtual loops of either SM particles or new states appearing in HV V interactions. The treatment of virtual photons has been illustrated with the JHU generator framework, and comparisons have been made to several other frameworks, including SMEFTsim using MadGraph5, HAWK, and Prophecy4f. A JHUGenLexicon program has been introduced for EFT parameterization and translation between frameworks. Overall, good agreement between the frameworks is found, including parameterization of EFT effects, once the sign conventions are matched, such as of the antisymmetric tensor ε µνρσ and the convariant derivative. The photon couplings appear to have a larger relative effect in decay compared to production, due to the dynamics of these processes. Nonetheless, the γH production topology appears interesting for isolating the photon couplings.
We have derived the scaling of the H boson production and decay rates with anomalous couplings, which are necessary for the total width calculations. From these, we also obtained the effective point-like couplings which reproduce the H → gg, γγ, and Zγ processes in the SM. We compared the effect of these couplings to the NLO EW corrections in the other processes with off-shell vector bosons. While the effects may be reproduced to within an order of magnitude, the effective couplings are not adequate for a careful modeling of the EW correction. At the same time, however, the squared higher-order terms may become as important as linear terms at lower q 2 values of the virtual photons, something that is neglected in the the NLO EW corrections. Subsequently, we make a proposal on how to handle singularities involving H boson decays to light fermions via photons. This procedure, based on matching amplitudes in the collinear limit, allows one to handle these singularities with efficient numerical evaluation. This approach can also incorporate the hadronic structure in decays to quarks using the experimentally measured quantities in low-energy processes.
We further make phenomenological observations on the special role of intermediate photons in analysis of LHC data in the EFT framework. Some of these features have been illustrated with projections for experimental measurements with the full LHC and HL-LHC datasets. The rates of the H → γγ and Zγ processes appear the most restrictive on the photon couplings, but cannot disentangle the CP -even and CP -odd couplings. We observe that the decay information H → 4 is more powerful in constraining the photon couplings than information in the VBF and V H production using the same H boson decay channel. However, these production channels, along with the γH production, can be analyzed using other H boson decays, with further gain in sensitivity to the photon couplings. The H → 4 channel can also be further optimized for the measurements of the photon couplings by relaxing the invariant mass and transverse momentum constraints on the leptons in analysis of LHC data. This will also require careful consideration of NLO EW effects in these processes, while including effects of anomalous contributions. Calculation of these effects will be required for precise analysis of HL-LHC data. add up to 1 due to interference effects. The ratio to the SM expectation is shown for C HX = 1. As discussed in more detail in Sections IV and VI, the presence of collinear singularities requires a special treatment of the low-q 2 V of the gauge bosons, and in the study in this Appendix we choose to apply a requirement q 2 V > 1 GeV 2 . The results of this study are discussed in Section III.      Table IV where q 2 ij,min = 4m 2 i , q 2 12,max = m 2 H and q 2 34,max = (m H − q 2 12 ) 2 . In this form, the two integrals over q 2 ij directly correspond to the V and V squared invariant masses of the matrix element M H→4f . This will be useful later.
In case of identical fermions, it is useful to symmetrize the phase space in the last line of Eq. (B5) because the matrix element has resonances not only in q 2 12 and q 2 34 but also in q 2 14 and q 2 32 , i.e. M = A(1234) + A(1432). Therefore, in practice we run the numerical simulation using dPS (4) where the integration boundaries are q 2 12,min = 4m 2 1 and q 2 12,max = m 2 H .

Collinear approximation applied to H → 2 γ
We split the invariant mass integration in Eq. (B8)