Jet azimuthal angle correlations in the production of a Higgs boson pair plus two jets at hadron colliders

Azimuthal angle correlations of two jets in the process $pp\to HHjj$ are studied. The loop induced $\mathcal{O}(\alpha_s^4 \alpha_{}^2)$ gluon fusion (GF) sub-process and the $\mathcal{O}(\alpha_{}^4)$ weak boson fusion (WBF) sub-process are considered. The GF sub-process exhibits strong correlations in the azimuthal angles $\phi_{1,2}^{}$ of the two jets measured from the production plane of the Higgs boson pair and the difference between these two angles $\phi_1^{}-\phi_2^{}$, and a very small correlation in the sum of them $\phi_1^{}+\phi_2^{}$. The impact of using a finite top mass $m_t^{}$ value on the correlations is found crucial. The transverse momentum of the Higgs boson can be used to enhance or suppress the correlations. The impact of a non-standard value for the triple Higgs self-coupling on the correlations is found much smaller than that on other observables, such as the invariant mass of the two Higgs bosons. The peak shifts of the azimuthal angle distributions reflect the magnitude of parity violation in the $gg\to HH$ amplitude and the dependence of the distributions on parity violating phases is analytically clarified. The WBF sub-process also produces correlated distributions and it is found that they are not induced by the quantum effect of the intermediate weak bosons but mainly by a kinematic effect. This kinematic effect is a characteristic feature of the WBF sub-process and is not observed in the GF sub-process. It is found that the correlations are different in the GF and in the WBF sub-processes. As part of the process dependent information, they will be helpful in the analyses of both the GF sub-process and the WBF sub-process at the LHC.


Introduction
The discovery of a Higgs boson with a mass around 125 GeV in 2012 is the main discovery of Run I of the Large Hadron Collider (LHC) [1,2]. The study of its properties has started and until now they are compatible with the Standard Model (SM) hypothesis [3][4][5]. To probe the mechanism of electroweak symmetry breaking (EWSB) [6][7][8][9] directly one would want to measure the triple Higgs self-coupling that is one of the key parameters of the scalar potential. This is one of the main goals of the future high-luminosity LHC and the Future Circular Collider (FCC) in hadron-hadron mode, a potential 100 TeV proton-proton collider following the LHC at CERN (for reviews of the FCC physics potential, see refs. [10,11]). In this view the production of a pair of Higgs bosons needs to be observed and has been extensively studied over the last years  (see also Refs. [33][34][35] for studies at the FCC). The gluon fusion (GF) sub-process [36][37][38][39] and the weak boson fusion (WBF) sub-process [36,[40][41][42] are the two main sub-processes, see e.g. Ref. [43] for a review of SM studies. Both of the two sub-processes are sensitive to the triple Higgs self-coupling. In the WBF sub-process, we have access to the coupling between the two Higgs bosons and the two weak bosons, too. In Refs. [23,31] studies using the production of a pair of Higgs bosons plus two hadronic jets has been conducted, using both the GF and the WBF sub-processes. The main advantage of the latter sub-process is the fact that the theoretical uncertainties are under control [19,28], but the phenomenological studies suffer from the difficulty to separate the GF contributions from that of the WBF.
Azimuthal angle correlations of two jets produced together with heavy particles have been actively studied as a provider of important information about the heavy particles [44][45][46][47][48][49][50]. The correlations are induced by only a certain type of sub-processes, called vector boson fusion (VBF) sub-processes, in which a heavy object is produced by a fusion of two vector bosons emitted from incoming two coloured particles. The correlations arise from the quantum effect of the two fusing intermediate vector bosons [42,47] 1 . A set of cuts on the rapidity y 1,2 of the two jets, y 2 < 0 < y 1 and y 1 − y 2 3 and an upper cut on the transverse momentum p T of the two jets are therefore crucial for the two jets to show strong correlations if any, since they enhance contributions from VBF sub-processes [47]. These rapidity cuts are often called VBF cuts.
The azimuthal angle correlations of the two jets both in the GF sub-process and in the WBF sub-process of the single Higgs boson plus two jets production process pp → Hjj [42,[44][45][46][47]52] are nowadays a common knowledge, have been studied in detail [53][54][55][56][57][58][59][60][61][62][63][64] and applied in many phenomenological studies. However, the azimuthal angle correlations of the two jets in the Higgs boson pair plus two jets production process pp → HHjj have not been studied thoroughly. To our knowledge, the correlations in the GF sub-process, which is an one-loop induced O(α 4 s α 2 ) process at leading order (LO), have not been studied in the literature. One of the reasons may be that event generations are still challenging even with an advanced calculation technique. For the GF sub-process in the process pp → Hjj, the approach of using the effective interactions between the Higgs boson and gluons is known to work quite well as long as the p T of the jets are small enough [65,66]. The azimuthal angle correlation after the VBF cuts can also be described correctly [66]. Therefore event generations can be easily performed with a tree-level event generator which implements the effective interactions. In contrast, for the GF sub-process in the process pp → HHjj, the effective interaction approach is known not to work well in describing the distributions of several observables [23,31]. It is naively expected that this observation is the same for the azimuthal angle correlations. The process pp → HHjj at LO with the exact one loop amplitude has been calculated for the first time in ref. [23] and subsequently phenomenology is studied in ref. [31]. The fully automated event generation for one-loop induced processes is now available in MadGraph5_aMC@NLO [67,68]. This achievement will activate further phenomenological studies of the process pp → HHjj including studies which use the azimuthal angle correlations. The azimuthal angle correlations of the two jets in the WBF sub-process, which is an O(α 4 ) process at LO, have been studied in [42]. There, only the azimuthal angle difference of the two jets is studied as an azimuthal angle observable.
In this paper, we study the azimuthal angle correlations of the two jets in the process pp → HHjj. Instead of considering all of the sub-processes contributing to the process pp → HHjj, we consider only the VBF sub-processes: gg → ggV * V * → ggHH (V = g), (1.1c) where q denotes a light quark or a light antiquark, g denotes a gluon and V * denotes an intermediate off-shell vector boson. The two Higgs bosons are produced by a fusion of two vector bosons emitted from incoming two coloured particles. More precisely, we calculate the amplitudes contributed from only t-channel Feynman diagrams shown in Figure 1 2 . We call them VBF amplitudes and VBF diagrams, respectively, in this paper. The circles denote all the Feynman diagrams contributing to the process V * V * → HH at LO. When the two virtual vector bosons are weak bosons (WBF sub-process), there are four tree-level Feynman diagrams, and when the two virtual vector bosons are gluons (GF sub-process), there are eight one-loop Feynman diagrams. We calculate these LO diagrams explicitly.
Our calculation is not exact but an approximated one. In ref. [47] it has been demonstrated that the inclusive cross section and the azimuthal angle correlations contributed only from the VBF diagrams are in good agreement with those contributed from all the diagrams, when the two jets are required to satisfy the VBF cuts and an upper p T cut, in the processes pp → Φjj, where Φ = H, A, G denotes the Higgs boson, a parity-odd Higgs boson and a spin-2 massive graviton, respectively. This observation indicates that a result based only on VBF diagrams is a good approximation to the exact LO result as long as the two jets satisfy the VBF cuts and an upper p T cut appropriately. Therefore our approach of calculating only the VBF amplitude is the simple and appropriate approach, when we study the azimuthal angle correlations (Let us remind that the correlations are induced only by the VBF sub-processes.). The same approach has been applied to the processes pp → QQjj, where Q denotes the top quark or the bottom quark, in ref. [50] and the validity of this approach is again observed. In order to check the validity of this approach in the GF sub-process by ourselves, we have numerically compared our cross section formula in which the loop-running top quark mass is set quite large (large m t limit) and the exact LO result generated by MadGraph5_aMC@NLO [67] which implements the effective interactions between gluons and the Higgs bosons (infinite m t limit). We have found a good agreement between the two results both in the inclusive cross section and in the distributions of several observables including azimuthal angle observables. From these observations, we are confident that our calculation of the azimuthal angle correlations is a good approximation to the exact LO result.
In order to measure the azimuthal angle correlations, we study four observables: φ 1 , φ 2 , ∆φ = φ 1 − φ 2 and φ + = φ 1 + φ 2 , where φ 1,2 are azimuthal angles of the two jets measured from the production plane of the Higgs boson pair (∆φ is irrelevant to the production plane as is clear from its definition). ∆φ is the observable which is sensitive to the property of the Higgs boson in the process pp → Hjj [44][45][46][47]49] and thus has been the subject of many studies. The processes pp → Gjj [47] and pp → QQjj [50] exhibit strong correlations in φ + . To our knowledge, correlations in φ 1,2 have not been addressed in any hadronic process in the literature. Our explicit findings can be summarised as follows. The GF sub-process has strong correlations in ∆φ and φ 1,2 , and the p T of the Higgs boson can be an useful measure to enhance or suppress these correlations. Using the finite m t value is important to produce the correlations correctly. Violation of the parity invariance of the gg → HH amplitude appears as the peak shifts of the correlations. The impact of a non-standard value for the triple Higgs self-coupling on the correlations is smaller than that on other observables, such as the invariant mass of the two Higgs bosons, of the inclusive process pp → HH. The correlation in φ + is negligibly small in most every case. The WBF sub-process produces correlated distributions in all of the azimuthal angle observables and they are not induced by the quantum effect of the intermediate weak bosons but mainly by a kinematic effect. This kinematic effect is a characteristic feature of the WBF sub-process and is not observed in the GF sub-process. The impact of a non-standard value for the triple Higgs self-coupling on the correlations is not significant in the WBF sub-process, too. The correlations in the GF and WBF sub-processes are found to be different.
The paper is organised as follows. In Section 2, we perform a calculation of the VBF amplitudes. Since it can be shown that the azimuthal angle correlations arise from the interference of amplitudes with various helicities of the two intermediate vector bosons [42,47], we employ a helicity amplitude technique. Our calculation is performed based on the method presented in ref. [47]. A full analytic set of helicity amplitudes is presented. Since the material in Section 2 is rather technical, the reader who is interested in only the results may skip this section. In Section 3, a detailed study of the azimuthal angle correlations is presented. First, we discuss the GF sub-process in Section 3.1. The squared VBF amplitude for the four sub-processes eq. (1.1) is given in a compact form. This analytic formula is found to be quite useful in making expectations of the correlations. The correlations in different kinematic regions of the two Higgs bosons and those in non-standard values for the Higgs triple self-coupling are studied. The impact of parity violation in the gg → HH amplitude on the correlations is also studied. Next, we discuss the WBF sub-process in Section 3.2. The squared VBF amplitude is given in a simple form by keeping only the dominant terms. The correlations in non-standard values for the Higgs triple self-coupling are studied. In Section 4 we summarise our findings and give some comments.
2 Helicity amplitudes for the process pp → V * V * jj → HHjj In this section we present a full analytic set of helicity amplitudes contributed from the vector boson fusion (VBF) diagrams. Our calculation is based on the method presented in ref. [47]. We present a more complete discussion on the treatment of the intermediate off-shell gluons in the VBF diagrams. In addition, we discuss the importance of an appropriate choice of gauge-fixing vectors for the polarisation vectors of the external gluons. We believe that the above two remarks provide sufficient justification for repeating some of the calculations of ref. [47].

VBF amplitudes
We first introduce a common set of kinematic variables where a 1,2,3,4 can be quarks, antiquarks or gluons, V 1,2 are intermediate off-shell vector bosons, H denotes the Higgs boson and the four-momentum k i and helicity σ i of each particle are shown. This assignment for the VBF sub-processes is more apparent in Figure 2. The external particles take helicities σ i = ±1 3 and the intermediate vector bosons take helicities λ i = ±1, 0. Colour indices are suppressed. The VBF helicity amplitude can be expressed as follows: 2) 3 We define the helicity operator for a two-component spinor by p · σ/| p| with the Pauli matrices σ, so a quark also takes σ = ±1 in our notation. Sometimes we simply write σ = ± instead of σ = ±1, and do the same for λ. Figure 2: The assignment of the four-momenta and helicities of each particle for the VBF sub-processes.
where J ) is a current involving the off-shell vector boson V i and the two external quarks, antiquarks or gluons, D is the tensor amplitude for the process V 1 V 2 → HH. When we denote a helicity amplitude in this paper, we show only helicity indices.
For the weak boson fusion (WBF) sub-process, we have (V 1 , V 2 ) = (W + , W − ), (W − , W + ) and (Z, Z), and only the VBF diagram shown in Figure 1 (left) contributes to the VBF amplitude. The VBF amplitude is gauge invariant on its own. This is apparent from the fact that only the VBF diagram exists in some of the WBF sub-processes, for instance there are no other diagrams describing the sub-process ud → udHH. We choose the unitary gauge for the weak boson propagator, We express the projector part of the above propagator in terms of polarisation vectors in the helicity basis (λ = ±1, 0, s): The λ = s component µ (q, λ = s) = q µ / −q 2 is the scalar part of a virtual weak boson and vanishes when it couples with a light quark current. As a result, the following replacement is possible when the external quarks are assumed to be massless: The polarisation vectors µ (q i , λ i = ±, 0) will be defined later once the kinematics of the weak bosons is fixed. After that, one can confirm eq. (2.4) explicitly. By inserting the identity eq. (2.5) into the VBF helicity amplitude in eq. (2.2), it can be expressed as a product of three helicity amplitudes: Eq. (2.7a) represents a helicity amplitude for the splitting process a i → a i+2 V i , where V i is off-shell. This will be derived in Section 2.2. Eq. (2.7b) represents a helicity amplitude for the process V 1 V 2 → HH, where V 1 and V 2 are off-shell. This will be presented in Section 2.3.
The gluon fusion (GF) sub-process (V 1 , V 2 ) = (g, g) is more complicated than the WBF subprocess. The VBF amplitude for the qq initiated sub-process (a 1 , a 2 ) = (q, q) is gauge invariant on its own, the reason being the same as for the WBF amplitude. If we choose the Feynman-'t Hooft gauge for a gluon propagator, the projector part of the propagator is: The λ = s component µ (q, λ = s) = q µ / −q 2 again vanishes when it couples with a quark current. As a result, the following replacement is possible without any approximation: Therefore, for the qq initiated sub-process, we can arrive at the same expression as in eq. (2.6). In ref. [47], the above procedure and thus the expression eq. (2.6) is used not only for the qq initiated sub-process but also for the qg and gg initiated sub-processes (a 1 , a 2 ) = (q, g), (g, g) 4 . We point out that this approach does not necessarily calculate the off-shell effects of the intermediate gluons correctly for the qg and gg initiated sub-processes and only introduces unnecessary complications in the amplitude calculation. Since the off-shell gluon amplitude (M HH gg ) λ 1 λ 2 is not enhanced in the on-shell limit of the gluons, we can always expand the off-shell gluon amplitude around the on-shell limit. To make our discussion simpler, let us consider an amplitude which involves only one off-shell gluon. The off-shell gluon amplitude can be expanded as where Q is the virtuality of the gluon and the first term in the right hand side (RHS) of the first equation is the on-shell amplitude, which is gauge invariant. As we will see in Section 2.2, the amplitude for the splitting process in eq. (2.7a), where an off-shell gluon with virtuality Q is emitted, has an overall factor of Q. By considering this factor and Q −2 in the propagator factor of the off-shell gluon, we find the following term in a VBF amplitude: (2.11) While the first term in the RHS of the above equation is gauge invariant, the rest terms are generally dependent on a gauge-fixing choice for the off-shell gluon. As we have mentioned above, for the qq initiated sub-process, the VBF amplitude is gauge invariant on its own and hence not only the first term in the RHS of eq. (2.11) but also the other terms as a whole are gauge invariant. This is not the case for the qg and gg initiated sub-processes. For the qg and gg initiated sub-processes, the second and higher terms in the RHS of eq. (2.11), which are not enhanced in the on-shell limit, become gauge invariant only after contributions from other diagrams are also included. Therefore, in our method where only the VBF diagrams are calculated, we can calculate the off-shell effect of the intermediate gluons correctly only for the qq initiated sub-process. For the qg and gg initiated sub-processes, therefore, we take the following approach: Completely ignore the off-shell effect of the intermediate gluons and look at only a kinematic region where the virtualities of the gluons are small (Q → 0). This is possible because the off-shell effect in the amplitude will not be essential as long as we look at the small virtuality region, as is clear from eq. (2.11). The unitarity condition gives the following equation for a gluon propagator in its on-shell limit q 2 → 0: By using this, we can express the VBF amplitude for the GF sub-process as a product of three helicity amplitudes in the same way that we do for the WBF sub-process and the qq initiated GF sub-process in eq. (2.6). The differences from the qq initiated GF sub-process are (1) ) λ 1 λ 2 is replaced by the on-shell gluon amplitude in which q 2 1 = 0 and q 2 2 = 0. Note that these two approximations are nothing less than the two fundamental ingredients of the equivalent photon approximation (EPA) [70,71]. As a rule of using the EPA, we should clarify the kinematic region where the approximated VBF amplitude is a good approximation to the exact VBF amplitude. For the process pp → HHjj, it should be −q 2 1 <ŝ and −q 2 2 <ŝ, whereŝ = M 2 HH . We use the EPA not only for the qg and gg initiated sub-processes but also for the qq initiated sub-process, so that we can consistently use the on-shell gluon amplitude (M HH . It should be noted that, while we use the on-shell gluon amplitude for the process gg → HH, the amplitude for an off-shell gluon emission in eq. (2.7a) is still calculated without ignoring the off-shell effect of the gluon (q 2 1,2 = 0), as in other applications of the EPA.

Helicity amplitudes for the splitting processes
We derive the helicity amplitude (2.7a) for the splitting processes. We use the chiral representation for Dirac matrices. Since we neglect the mass of the external quarks, the helicity of the quark is equal to its chirality and the helicity of the antiquark is opposite to its chirality. As a result, the helicity of a i is always equal to that of a i+2 (σ i = σ i+2 ) for the quark and antiquark splitting processes. Otherwise, the amplitude vanishes. By introducing one common currentĴ µ i , we can express the quark and antiquark currents J µ where u(k, σ) α represents the two-component Weyl u-spinor with its four-momentum k, helicity σ and chirality α (= ±1), and σ µ ± = (1, ± σ) with the Pauli matrices σ. Note that we can use the aboveĴ µ i for the antiquark current, too, since (2.14) The couplings between quarks and vector bosons relevant to our study are summarised as follows: where g s is the QCD coupling, t a is the generator of the SU(3) group, g w = e/ sin θ w , g z = e/(sin θ w cos θ w ) with θ w being the electroweak mixing angle and e being the proton charge, V ud is an element of the CKM matrix, Q q is the electric charge of the quark in unit of e, T 3 u = 1/2 and T 3 d = −1/2. The gluon current involving three gluons is also expressed by eq. (2.13a) with where f abc is the structure constant of the SU(3) group.
Generally speaking, helicities are frame dependent and so are helicity amplitudes. When we calculate the VBF helicity amplitude M σ 3 σ 4 σ 1 σ 2 in eq. (2.6), we must choose one frame at first and then define the four-momenta and the helicities of all the external particles in this particular frame. For our calculation we choose the centre-of-mass (c.m.) frame of the two intermediate vector bosons moving along the z-axis, which is shown in the middle of Figure 3. We call it the VBF frame. Note that the production plane of the two Higgs bosons coincides with the plane of the x-z axes. All of the three helicity amplitudes in the VBF helicity amplitude must be evaluated in the VBF frame. However, by using a property of helicities, we can justify a calculation of the helicity amplitude (J for the splitting processes in a different frame. Because the helicity of a massless quark is frame independent and furthermore the helicity of a massive vector boson is invariant under Lorentz boosts along its momentum direction as long as the boosts do not change the sign of its three-momentum 5 , the helicity amplitude (J for the quark splitting processes is invariant under Lorentz boosts along the z-axis from the VBF frame, as long as the boosts do not change the sign of the three-momentum of V i . By this property, it is justified to calculate the helicity amplitude (J for the quark splitting process in the q 1 Breit frame, to which k 1,3 and q 1 in the VBF frame can move by a single Lorentz boost along the negative direction of the z-axis. Similarly, it is justified to calculate the helicity amplitude (J V 2 a 2 a 4 ) λ 2 σ 2 σ 4 for the quark splitting process in the q 2 Breit frame, to which k 2,4 and q 2 in the VBF frame can move by a single Lorentz boost along the positive direction of the z-axis. The q 1 and q 2 Breit frames are illustrated in the left and right of Figure 3, respectively. A calculation of the helicity amplitude (J for the gluon splitting process in the q i Breit frame is also justified by appropriately choosing gauge-fixing vectors for the polarisation vectors of the external gluons. Although the helicity of a gluon is frame independent, the polarisation vectors of the gluon are dependent on their gauge-fixing vectors and hence the helicity amplitude for the gluon splitting process is in general frame dependent. However, if we choose the gauge fixing vectors in a way that their directions are invariant under Lorentz boosts along the z-axis, the helicity amplitude for the gluon splitting process also becomes invariant under the same boosts, again as long as the boosts do not change the sign of the three-momentum of the intermediate off-shell gluon. This can be simply achieved by choosing all the gauge-fixing vectors along the z-axis.
We parametrise the four-momenta k 1 , k 3 and q 1 in the q 1 Breit frame as where Q 1 = −(q 1 ) 2 > 0, 0 < θ 1 < π/2 and 0 < φ 1 < 2π, and the four-momenta k 2 , k 4 and q 2 in the q 2 Breit frame as where Q 2 = −(q 2 ) 2 > 0, π/2 < θ 2 < π and 0 < φ 2 < 2π. We define polarisation vectors for the intermediate vector boson V 1 in the q 1 Breit frame by and those for the intermediate vector boson V 2 in the q 2 Breit frame by It is easy to explicitly confirm that the above sets of the polarisation vectors satisfy the identity for the propagator of a weak boson in eq. (2.4). We use the same polarisation vectors µ (q i , λ i = ±) for the intermediate weak bosons for a weak boson emission from a quark and those for a gluon emission from a quark are the same, except for couplings. (Let us remind that the λ = 0 components of the intermediate gluons are neglected.) As we have discussed above, the gauge-fixing vectors for polarisation vectors of the external gluons should be chosen along the z-axis. For the polarisation vectors of the two external gluons in the q 1 Breit frame, we choose a common vector n µ 1 = (1, 0, 0, −1), i.e. the light-cone axial gauge. With this choice, the polarisation vectors in the q 1 Breit frame are Similarly, a vector n µ 2 = (1, 0, 0, 1) is commonly chosen for the polarisation vectors of the two external gluons in the q 2 Breit frame, then the polarisation vectors in the q 2 Breit frame are It is easy to confirm that the above polarisation vectors with their gauge fixing vectors satisfy the unitarity condition for an on-shell gluon: With our preparations up to now, we can easily derive the helicity amplitude (J for the quark splitting process and that for the gluon splitting process in the Breit frames. Following ref. [47], we write the amplitudes as (2.24)  The coupling factor is already defined in eq. (2.13a). The common amplitudeĴ is summarised in Table 1. Note that we adopt the phase convention for the two-component Weyl spinors developed in refs. [72,73]. Since we use the same phase convention for the spinors and the same gauge fixing vectors for the external gluons with ref. [47], the amplitudes in Table 1 are consistent with those in tables 1 and 2 of ref. [47] including the common overall phase.

Helicity amplitudes for the processes V V → HH
Finally we present the helicity amplitudes for the processes defined in eq. (2.7b). As we have mentioned in Section 2.2, we evaluate the amplitudes in the VBF frame. We parametrise the four-momenta q 1 , q 2 , q 3 and q 4 in the VBF frame as whereŝ is the c.m. energy squaredŝ = (q 1 + q 2 ) 2 and β = 1 − 4m 2 H /ŝ with m H being the Higgs boson mass. The polarisation vectors for the vector bosons V 1,2 in the VBF frame can be simply obtained by boosting those in the q 1,2 Breit frames (eqs. (2.19) and (2.20)) to the VBF frame along the z-axis: and Needless to say, the λ 1,2 = ±1 components remain the same after the boost. The helicity amplitude for the WBF process , is given by [42] M HH where G F is the Fermi constant and λ h is the factor re-scaling the triple Higgs self-coupling: The standard model predicts λ h = 1. The notation for the propagators is defined in eq. (2.3). By using the four-momenta and the polarisation vectors defined above, we can obtain the off-shell weak boson amplitude.
As we have discussed in Section 2.1, we use the on-shell gluon amplitude for the GF process gg → HH. Using the notation given in the appendix of ref. [39], we write the amplitude as where b 1,2 are the colour indices of the gluons, and F , F and G are form factors [39] consisting of the scalar loop functions after tensor reduction and A iµ 1 µ 2 are the tensor structures of the process. Not only we neglect the λ 1,2 = 0 components of the gluons, but we also set Q 1 = Q 2 = 0 in the four-momenta in eq. (2.25). Our definitions of the polarisation vectors µ (q i , λ i = ±) in eqs. (2.26) and (2.27) actually correspond to an axial gauge n µ = (1, 0, 0, −1) for an on-shell gluon q µ 1 ∝ (1, 0, 0, 1) and to an axial gauge n µ = (1, 0, 0, 1) for an on-shell gluon q µ 2 ∝ (1, 0, 0, −1), respectively.
In order to simplify further analyses, we introduce the following amplitudê Then, the amplitude eq. (2.29) has a simpler form After a simple manipulation in the VBF frame, we find The parity invariance of the amplitude is apparent, (M HH gg ) ++ = (M HH gg ) −− and (M HH gg ) +− = (M HH gg ) −+ . Up to now we have assumed the Higgs sector of the standard model. In this paper, we study the impact of parity symmetry violation in the GF process gg → HH 6 . We can read off parity violating phases from a parity-odd process gg → HA, where A is a parity-odd Higgs boson. By using the tensor A iµ 1 µ 2 for the process gg → HA [39], we find the amplitude in the VBF frame: In order to evaluate the impact of parity violation, we introduce two angles (or phases) ξ 1,2 (−π/2 ≤ ξ 1,2 ≤ π/2) and write the amplitude aŝ The two phases ξ 1,2 parametrise the magnitude of parity violation in the process gg → HH and independently affect the ∆λ = λ 1 − λ 2 = 0 helicity states and ∆λ = ±2 helicity states, respectively. We believe that this simplified introduction of parity violating phases is enough for studying the impact of parity violation on the azimuthal angle correlations.

Azimuthal angle correlations
In this section we present a detailed study of the azimuthal angle correlations of the two jets by using the helicity amplitudes presented in Section 2.

The gluon fusion process
The azimuthal angle correlations of the two jets can be analytically apparent, once we obtain the squared VBF amplitude. There are four gluon fusion (GF) sub-processes (V 1 , V 2 ) = (g, g), namely the qq initiated sub-process (a 1 , a 2 ) = (q, q), the qg initiated sub-processes (a 1 , a 2 ) = (q, g), (g, q) and the gg initiated sub-process (a 1 , a 2 ) = (g, g). The squared VBF amplitude for the four GF 6 The parity symmetry violation in the GF process indicates charge-conjugation and parity (CP) symmetry violation in the Higgs sector.  sub-processes has the following compact form, after averaging over the initial state colours and helicities and summing over the final state colours and helicities: In order to simplify our writing, we introduce a notationM λ 1 λ 2 which denotes the helicity amplitude (M HH gg ) λ 1 λ 2 . C a 1 a 2 are the colour factors from the splitting processes a 1 → a 3 + g * 1 and a 2 → a 4 + g * 2 , and they take values C qq = 16/9, C qg = C gq = 4 and C gg = 9. The colour factors C a 1 a 2 and the functions F i [a 1 a 2 ] for antiquarks are the same as those for quarks, since the QCD interaction does not distinguish quarks and antiquarks.
The first term in the right hand side (RHS) of eq. (3.1) contributes to the inclusive cross section after a phase space integration, while the other terms give the azimuthal angle distributions of the two jets. The azimuthal angles φ 1,2 of the two jets are defined in the q 1,2 Breit frames, respectively, and they remain the same in the VBF frame. In the limit that each of the two jets in the protonproton (pp) frame is collinear to the incoming parton that emits it (collinear limit), the emitted two intermediate vector bosons also move on the z-axis. After rotating the two jets and the two Higgs bosons around the z-axis in such a way that the two Higgs bosons have zero azimuthal angle (Let us remind that the two Higgs bosons have zero azimuthal angle in the VBF frame, see eq. (2.25)) and after a single boost along the z-axis, all of these particles can be studied in the VBF frame. Therefore, in the collinear limit, the azimuthal angles of the two jets after the single rotation around the z-axis are identical to φ 1 and φ 2 . We apply the VBF cuts and an upper transverse momentum p T cut on the jets in the pp frame and these cuts reproduce the collinear limit to some extent. Hence the azimuthal angles of the two jets in the pp frame after the rotation around the z-axis should not be very different from φ 1,2 defined in the q 1,2 Breit frames. We perform the rotation of the two jets and the two Higgs bosons around the z-axis in the following way: 1. Go to the centre-of-mass (c.m.) frame of the two Higgs bosons and then rotate the two Higgs bosons around the z-axis byφ in a way that the two Higgs bosons have zero azimuthal angle.
2. Rotate the two jets in the pp frame around the z-axis byφ.
3. Measure the azimuthal angles of the two jets.
In the collinear limit, it is clear that the azimuthal angles of the two jets measured after this rotation coincide with φ 1,2 . Note that this rotation is necessary for the azimuthal angles and the sum of them to show meaningful distributions, because the process pp → HHjj in the pp frame is completely symmetric around the z-axis. This rotation is, however, not needed for the difference of the two azimuthal angles. Before we show numerical results, we point out characteristic features of the GF sub-process in the standard model (SM) already expected from the analytic formula eq. (3.1): • The azimuthal angles of the two jets show the same distribution due to the parity invariance of the amplitudeM λ 1 λ 2 , see eq. (2.32).
• All of the azimuthal angle observables show cosine distributions, again due to the parity invariance of the amplitude.
Violation of the parity invariance of the amplitude should appear as a deviation from the above expectations. This case will be studied at the end of this subsection.
We show numerical results for the 14 TeV LHC. We do not study decays of the two Higgs bosons and assume that they can be reconstructed. An outgoing quark, antiquark or gluon is identified as a jet. The following set of parameters are chosen: m H = 125.5 GeV, m t = 173.5 GeV and α s (m Z ) = 0.13. We use the CTEQ6L1 [74] set for the parton distribution functions (PDFs) and have chosen the input value for the strong coupling constant accordingly. For the scales in the PDFs, we choose a fixed value of 25 GeV, which corresponds to the lower cutoff on the transverse momentum p T of the jets (see below). The scales in the strong couplings are chosen as α s ( √ŝ ) 2 α s (150GeV) 2 , where √ŝ is the invariant mass of the two Higgs bosons. Using the two different scales in the strong couplings can be considered as a better choice, since we look at only a kinematic region where the virtualities of the gluons are small (Q 1,2 → 0) and this separates the two splitting processes from the process gg → HH in time-scale. The scales in the strong couplings of the splitting processes correspond to the upper cutoff on the p T of the jets. The following cuts are applied on the rapidity y and p T of the two jets in the pp frame: The above rapidity cuts are the VBF cuts. As we have already mentioned in Section 1, the VBF cuts and the upper p T cut enhance the contributions from the VBF sub-processes. This can be understood as follows. The virtuality Q 1 of the intermediate vector boson V 1 is where the momentum assignment given in eq. (2.1) is used and θ 13 is the angle between k 1 and k 3 . The VBF cuts make k 3 collinear to k 1 (θ 13 → 0), then Q 1 is decreased and the VBF amplitude is enhanced. The upper p T cut additionally enhances the VBF amplitude. In a collinear case (θ 13 → 0), the p T of k 3 is so the upper p T cut reasonably implies the upper cut on Q 1 . Only the VBF cuts may be enough to enhance contributions from the VBF sub-processes. However, we need to impose the upper p T cut, too, because we perform the two approximations in calculating the VBF amplitude for the GF sub-process and these approximations are justified only when Q 2 1 <ŝ and Q 2 2 <ŝ. Throughout our analyses, the four-momentum of the jet which has a positive rapidity y 1 in the pp frame is used for calculating its azimuthal angle labelled as φ 1 and that of the other jet which has a negative rapidity y 2 in the pp frame is used for calculating its azimuthal angle labelled as φ 2 . Azimuthal angles labelled as φ 1,2 in our numerical results shown below are not those defined in the q 1,2 Breit frames anymore. The phase space integration and event generations are performed with the programs BASES and SPRING [75]. The scalar loop functions are calculated with the program FF [76].
In Figure 4 we show the normalised differential cross section as a function of φ 1 (upper left), φ 2 (upper right), ∆φ = φ 1 − φ 2 (lower left) and φ + = φ 1 + φ 2 (lower right). The blue solid curve, labelled as (1) VBF, represents the result according to our analytic cross section formula, to which only the VBF diagrams contribute. The scalar loop functions in the form factors are calculated by using the finite m t value. The red solid curve, labelled as by implementing the following effective Lagrangian density [77] into an UFO file [78] where F a µν is the gluon field strength tensor and H is the Higgs boson field. The interactions between gluons and the Higgs bosons in this approach can be considered as those induced by a quark loop where the mass of the quark is infinitely large. The consistency between the result (2) and the result (3) in all of the panels shows that our analytic cross section formula is a good approximation to the exact cross section formula. The discrepancy between the result (1) and the result (2) in the φ 1,2 plots particularly shows the importance of using a finite m t value in the azimuthal angle distributions. The reason why we do not find correlated distributions in φ 1,2 in the large m t results can be understood from the squared VBF amplitude in eq. (3.1). The second and third terms in the RHS shows that the correlations in φ 1,2 arise from the interference  of the amplitudesM λ 1 λ 2 for ∆λ = 0 and ∆λ = ±2 states, where ∆λ = λ 1 − λ 2 . In the large m t limit, the amplitude for ∆λ = ±2 states vanishes and hence the correlations also vanish. The GF sub-process exhibits the largest correlation in ∆φ and an almost zero correlation in φ + . The squared VBF amplitude in eq. (3.1) tells us that the correlation in ∆φ arises from the interference of the amplitudesM λ 1 λ 2 for ∆λ = 0 states (the fourth term in the RHS) and the correlation in φ + arises from the interference of the amplitudes for ∆λ = ±2 states (the fifth term in the RHS). The reason of the large correlation in ∆φ and the small correlation in φ + is because the amplitudes for ∆λ = 0 states are much larger than those for ∆λ = ±2 states in a large part of the phase space (this will be confirmed explicitly below soon). We briefly mention differences between the processes pp → HHjj and pp → Hjj in the correlations. The process gg → H has non-zero amplitudesM λ 1 λ 2 only for ∆λ = 0 states. Therefore, the process pp → Hjj exhibits a large correlation in ∆φ [44][45][46][47] but no correlations in φ 1,2 and φ + .
In Figure 5 we show the differential cross section as a function of φ 1,2 (left), ∆φ = φ 1 − φ 2 (middle) and φ + = φ 1 + φ 2 (right), contributed by all of the sub-processes (black curve, labelled as (1) pp), by the qg and gq initiated sub-processes (blue curve, labelled as (2) qg + gq), by the gg initiated sub-process (red curve, labelled as (3) gg) and by the qq initiated sub-process (green curve, labelled as (4) qq). All of the numerical results hereafter in this subsection are produced by using our analytic cross section formula. As we have already discussed before showing the numerical results and have actually confirmed in Figure 4, φ 1 and φ 2 show the same distribution due to the parity invariance of the amplitude in the SM. Therefore, we show only the φ 1 distribution and label it φ 1,2 instead of showing the both distributions, until we study parity violation. In Table 3, we present the inclusive cross sections of qqHH, ggHH, qgHH final states and the sum of these final states in pp collisions. The calculation methods are the same as those used in Figure 4. A good agreement (within 5%) between the second row and the third row for each cross section again confirms the validity of our analytic cross section formula. A large discrepancy between the first row and the second (or third) row for each cross section can be considered as another important result of using the finite m t value.
In order to study how the azimuthal angle correlations depend on the kinematics of the two  Higgs bosons, we introduce the following three quantities: These are simply the coefficients of the azimuthal angle dependent terms divided by the coefficient of the azimuthal angle independent term in eq. (3.1). Although the functions F i [a 1 a 2 ] are omitted, these quantities are useful in that an azimuthal angle observable should show a strong correlation in a kinematic region where the corresponding quantity CR i is enhanced. We find that all of the quantities CR i are well correlated with the p T of the Higgs boson in the c.m. frame of the two Higgs bosons, as shown in Figure 6. Figure 6 presents list plots of CR 2φ 1,2 (left), CR 2∆φ (middle) and CR 2φ + (right) with the p T of the Higgs boson in the c.m. frame of the two Higgs bosons. The large value of CR 2∆φ and the small value of CR 2φ + in the most part of the phase space explain the large correlation in ∆φ and the small correlation in φ + observed in Figure 4. The plots indicate that the correlation in ∆φ is enhanced as the p T of the Higgs boson is decreased, while the correlations in φ 1,2 and φ + are enhanced as the p T of the Higgs boson is increased. This is confirmed in Figure 7. In Figure 7, we show the normalised differential cross section as a function of φ 1,2 (left), ∆φ = φ 1 − φ 2 (middle) and φ + = φ 1 + φ 2 (right) with different values for the lower cutoff on the p T of the Higgs boson in the c.m. frame of the two Higgs bosons, blue solid curve: p T > 0 GeV (no cut), black dashed curve: p T > 150 GeV, and red solid curve: p T > 250 GeV.
Next, we study how the azimuthal angle correlations depend on the triple Higgs self-coupling. Eqs. (2.31) and (2.32) tell us that λ h , which is the factor re-scaling the triple Higgs self-coupling, affects only the amplitudeM λ 1 λ 2 for ∆λ = 0 states. However, as is clear from the quantities For an extreme example, if we could make the amplitude for ∆λ = 0 states to be identically zero (this is actually not possible because the amplitude depend not only on λ h but also on the kinematic of the process gg → HH), we would observe a large correlation only in φ + . In Figure 8, we study the correlations with three different values for λ h , the blue solid curve: λ h = 0, the black dashed curve: λ h = 1 (the SM prediction), and the red solid curve: λ h = 2. The impact of a non-standard value for λ h (λ h = 1) in the distributions is visible but not large. Actually this is much smaller than that in other observables, such as the p T of the Higgs boson [18,19] or the invariant mass of the two Higgs bosons [19,22], of the inclusive process pp → HH. The azimuthal angle correlations may not be useful to probe λ h .
Finally, we study the impact of parity symmetry violation on the azimuthal angle correlations. Let us remind that we have introduced two phases ξ 1,2 in eq. (2.34) and these phases parametrise the magnitude of parity violation in the process gg → HH. If they are non-zero, the gg → HH amplitudeM λ 1 λ 2 is not parity invariant anymore:M ++ =M −− if ξ 1 = 0,M +− =M −+ if ξ 2 = 0. The squared VBF amplitude eq. (3.1) tells us that violation of the parity invariance of the amplitude appears as the following deviations from the standard model predictions: (1) the φ 1 and φ 2 distributions are not necessarily equal to each other, (2) the azimuthal angle observables do not necessarily show cosine distributions. In order to make these expectations more explicit, we write the amplitudeM λ 1 λ 2 with the phases ξ 1,2 in the following way: M λ,−λ = B e +iλξ 2 .
These results show that parity violation appears as peak shifts of the azimuthal angle distributions. The results also tell us an important fact that the peak shifts of the distributions reflect only the magnitude of parity violation (Recall that ξ 1,2 parametrise the magnitude of parity violation in our study). It is an easy exercise to confirm that this is true no matter how parity violating phases are introduced. From the above analytic results, it is also apparent how each observable depends on the phases ξ 1,2 . ∆φ is sensitive to ξ 1 and φ + is sensitive to ξ 2 . φ 1 and φ 2 are sensitive to both of ξ 1 and ξ 2 . Although we are far less likely to be able to measure ξ 2 in the φ + distribution because of the very small correlation in φ + , we may use φ 1 and φ 2 to probe ξ 2 instead.
We note in passing the impact of parity violation in the process pp → Hjj, which can be considered as a simpler case. Since the process gg → H has non-zero amplitudesM λ 1 λ 2 only for ∆λ = λ 1 − λ 2 = 0 states, the process pp → Hjj has only eq. (3.9c) in its cross section formula. Therefore, parity violation in the process gg → H appears as a peak shift only in the ∆φ distribution [45,46].
While the phase dependence on the correlations is apparent from eq. (3.9), we present numerical results, too. In Figure 9 we show the normalised differential cross section as a function of φ 1 (left column), φ 2 (middle column) and ∆φ = φ 1 − φ 2 (right column) with different values for ξ 1,2 . The correspondence between the curves and values for ξ 1,2 is shown inside each panel. In all of the panels, the SM prediction is shown by the black dashed curve. A cutoff, p T > 200 GeV, is imposed on the p T of the Higgs boson in the c.m. frame of the two Higgs bosons, when we produce the φ 1 and φ 2 plots, in order to enhance the correlations in φ 1,2 , see Figure 7. The distributions of φ + = φ 1 + φ 2 are not shown anymore, since we have found that the correlation in φ + is very small in most every case, see Figures 7 and 8.

The weak boson fusion process
The weak boson fusion (WBF) sub-process (V 1 , V 2 ) = (W + , W − ), (W − , W + ) and (Z, Z) consists of only the qq initiated sub-process (a 1 , a 2 ) = (q, q). However, the squared VBF amplitude for the WBF sub-process takes a more complicated form than that for the GF sub-process in eq. (3.1) because of the following two reasons: (1) the helicity λ 1,2 = 0 components of the intermediate weak bosons additionally induce eight azimuthal angle dependent terms, such as cos (φ 1 − φ 2 ), (2) we cannot simply take averages for the initial helicities and take summations for the final helicities, since the electroweak interactions distinguish different helicity states. The squared VBF amplitude for four sets of the helicities of the external quarks must be prepared: The amplitude (M HH for helicity λ 1,2 = 0 weak bosons gives a dominant contribution in the WBF sub-process and the amplitudes for other helicities are much smaller than the above amplitude. Thus we can expect that all of the azimuthal angle correlations are small, because the correlations arise from the interference of the amplitudes for various helicities, see eq. (3.1). If we keep only terms which contain at least one (M HH ) 00 , the squared VBF amplitude for (σ 1 , σ 3 , σ 2 , σ 4 ) = (+, +, +, +) has the following form: where we introduce a notationM λ 1 λ 2 which denotes the amplitude (M HH The squared VBF amplitude for the other three helicity states can be simply obtained by exchanging s 1,2 and t 1,2 in |M ++ ++ | 2 in the following ways: The couplings should be also changed accordingly. The coefficients of cos φ 1,2 terms actually contain oneM 00 , too. However, the cos φ 1,2 terms cannot give correlated distributions in the process pp → HHjj, because we cannot distinguish the two Higgs bosons. More practically, an azimuthal angle dependent term which changes its overall sign under the transformations φ 1 → φ 1 + π and φ 2 → φ 2 + π gives only a flat distribution after the phase space integration. The first term in the right hand side (RHS) of eq. (3.10) contributes to the inclusive cross section after the phase space integration and the other terms give the correlations in ∆φ = φ 1 − φ 2 and φ + = φ 1 − φ 2 . An interesting difference from the GF sub-process is that the sine terms do not vanish even when the amplitude is parity invariant (M ++ =M −− andM +− =M −+ ), if the amplitude contains an imaginary part. This is because the interaction between the external quarks and the intermediate weak boson already violates the parity symmetry. In our tree-level calculation, the amplitude is purely real and we will observe only cosine distributions.
We show numerical results for the 14 TeV LHC. The setup and phase space cuts are the same as in the GF study in Section 3.1. Therefore the numerical results in this Section can be directly compared with those in Section 3.1. The only difference is the scale choice in the PDFs. The p T of the jet with a positive rapidity is used for the scale in the PDF of the incoming parton which moves along the positive direction of the z-axis, and the p T of the jet with a negative rapidity is used for the scale in the PDF of the other incoming parton. In Figure 10, we show the normalised differential cross section as a function of φ 1 (upper left), φ 2 (upper right), ∆φ = φ 1 − φ 2 (lower left) and φ + = φ 1 + φ 2 (lower right). The blue solid curve, labelled as (1) (3) in the ∆φ and φ + plots come from the cos ∆φ and the cos φ + terms in eq. (3.10), respectively. The smallness of the discrepancies is as expected (see the discussion above eq. (3.10)). In the ∆φ and φ + plots, the result (3) again shows correlated distributions, which must be induced by a kinematic effect. Therefore, the WBF sub-process produces the correlated distributions in all of the azimuthal angle observables, however they are mainly induced by a kinematic effect. We note that the GF sub-process produces only flat distributions in all of the azimuthal angle observables, when the quantum effects of the intermediate gluons are killed in the same way as above. It can be concluded that the non-flat distributions induced by a kinematic effect are a characteristic feature of the WBF sub-process.
We study how the azimuthal angle distributions depend on the triple Higgs self-coupling. We have observed the kinematic effect on the distributions in non-standard cases λ h = 1, too, where λ h is the factor re-scaling the triple Higgs self-coupling. The three cases λ h = 0, 1, 2 produce the similar distributions in the azimuthal angle observables, when the quantum effects of the intermediate weak bosons are killed. In Figure 11 we show the normalised differential cross section as a function of φ 1,2 (left), ∆φ = φ 1 − φ 2 (middle) and φ + = φ 1 + φ 2 (right) with three different values of λ h , the blue solid curve: λ h = 0, the black dashed curve: λ h = 1 (the SM prediction), and the red solid curve: λ h = 2. The correlated distributions in the φ 1,2 plot are completely induced by a kinematic effect. We find that, when λ h = 2, the coefficient of the cos (φ 1 − φ 2 ) term (the second term in the RHS of eq. (3.10)) is large enough to flip the ∆φ distribution. The impact of a non-standard value for λ h (λ h = 1) is again not so significant. However, differently from the GF sub-process which is actually the O(α 2 s ) correction to the inclusive GF sub-process gg → HH, the WBF sub-process is a LO tree-level process and so the correlations should be used together with other observables, such as the invariant mass of the Higgs boson pair, to probe λ h .   Figure 11: The normalised differential cross section of the WBF process as a function of φ 1,2 (left), ∆φ (middle) and φ + (right), with three different values for the triple Higgs self-coupling re-scaling factor λ h .
The correspondence between the curves and values for λ h is shown inside the left panel.

Summary and discussion
In this paper, we have studied the azimuthal angle correlations of two jets in the production of a Higgs boson pair plus two jets pp → HHjj. Based on the known fact that the azimuthal angle correlations are induced by the quantum effects of the two intermediate vector bosons in vector boson fusion (VBF) sub-processes, we have calculated the amplitudes contributed from only VBF Feynman diagrams. As VBF sub-processes, we have considered the gluon fusion (GF) sub-process, which is an one-loop O(α 4 s α 2 ) process at leading order (LO), and the weak boson fusion (WBF) sub-process, which is an O(α 4 ) process at LO. We have used a helicity amplitude technique for evaluating the VBF amplitudes. Based on the method presented in ref. [47], we have divided a VBF amplitude into two amplitudes for off-shell vector boson emissions (q → qV * or g → gV * ) and one amplitude for the Higgs boson pair production by a fusion of the two off-shell vector bosons (V * V * → HH), and presentd each of the three amplitudes in the helicity basis. The quantum effects of the intermediate vector bosons are still included correctly and are expressed as the interference of the V * V * → HH amplitudes with various helicities of the vector bosons. With this method, we have obtained the analytic cross section formula in a compact form, from which we can easily make an expectation on the azimuthal angle correlations, both for the GF sub-process and for the WBF sub-process. We have numerically compared our analytic cross section formulas with the exact LO results and have observed the good agreement between the two results both in the inclusive cross sections and in azimuthal angle distributions, after the VBF cuts and the upper transverse momentum p T cut on the jets are imposed (For the GF sub-process, the comparison is performed only in the large m t limit.). As azimuthal angle observables, we have studied four observables: φ 1 , φ 2 , ∆φ = φ 1 − φ 2 and φ + = φ 1 + φ 2 , where φ 1,2 are the azimuthal angles of the two jets measured from the production plane of the Higgs boson pair.
In the GF sub-process, using a finite m t value is found to be important to produce the azimuthal angle correlations correctly. The GF sub-process exhibits large correlations in φ 1,2 and ∆φ. The p T of the Higgs boson is found to be useful in controlling these correlations. The correlation in ∆φ is enhanced when the p T of the Higgs boson is decreased and the correlations in φ 1,2 are enhanced when the p T of the Higgs boson is increased. We have found that the correlation for the GF process (red solid curve) and the WBF process (black dashed curve).
in φ + is very small in most every case. The impact of a non-standard value for the triple Higgs self-coupling on the correlations is found to be much smaller than that in other observables, such as the invariant mass of the Higgs boson pair, of the inclusive process pp → HH. In order to study the impact of parity violation on the correlations, we have introduced two independent phases ξ 1,2 which parametrise the magnitude of parity violation in the process gg → HH, in a way that ξ 1 affects the gg → HH amplitude for ∆λ = 0 helicity states, where ∆λ = λ 1 − λ 2 and λ 1,2 are helicities of the gluons, and ξ 2 affects the amplitude for ∆λ = ±2 helicity states.
We have analytically shown that parity violation appears as peak shifts of the correlations and that the peak shifts reflect only the magnitude of parity violation. We have also shown that ∆φ is sensitive to ξ 1 , φ + is sensitive to ξ 2 , and φ 1,2 are sensitive to both of ξ 1 and ξ 2 . Although we are far less likely to be able to measure ξ 2 in the φ + distribution because of the very small correlation in φ + , we may use φ 1 and φ 2 to probe ξ 2 instead. While we can naively expect that the azimuthal angle distributions in the WBF sub-process are almost flat, we have actually observed correlated (non-flat) distributions. We have found that they are not induced by the quantum effects of the intermediate weak bosons but by a kinematic effect. Since we do not find the similar kinematic effect in the GF sub-process, we conclude that this is a characteristic feature of the WBF sub-process. The impact of a non-standard value for the triple Higgs self-coupling is not so significant in the WBF sub-process, too.
The parton level event samples of the process pp → HHjj are exclusively generated and each of the two outgoing partons is identified as a jet in our numerical studies. When more realistic event generations are intended to be performed, merging the parton level event samples with the leading logarithmic parton shower [80][81][82] and subsequently proceeding to a hadronisation procedure will be a promising approach. However, a careful merging procedure is required for correctly reproducing the azimuthal angle correlations after the merging procedure, because the correlations studied in this paper are completely process dependent (pp → HHjj) and those process dependent angular correlations are not described correctly by the parton shower. Contamination from the parton shower may lead to a wrong prediction [83].
The azimuthal angle correlations revealed in this paper will help the analyses of the process pp → HHjj at the LHC. Since the production cross section of the process pp → HHjj is small, we should use as much process dependent information as possible to extract the events of the process. The azimuthal angle correlations are obviously a part of the process dependent information. It is a known issue that separating the contributions coming from the GF sub-process and those coming from the WBF sub-process is difficult in the production of a Higgs boson pair plus two jets [23]. One possible application of the correlations is to help disentangling the GF sub-process and the WBF sub-process by using the fact that these two sub-processes exhibit different correlations, as shown in Figure 12. We have studied only the signal processes in this paper and so the impact of the correlations in a realistic situation is not so clear yet. The fully automated event generation for loop induced processes is now available in MadGraph5_aMC@NLO [67,68]. This achievement will activate phenomenological studies of the process. We hope that further phenomenological studies including the uses of the azimuthal angle correlations will be performed both by theorists and by experimentalists.