$e^{+}e^-\rightarrow hhZ$ in the B-L symmetric SSM

The double Higgs boson production through $e^{+}e^-\rightarrow hhZ$ is analyzed in the minimal supersymmetric extension of the Standard Model(SM) with the local gauge symmetry $U(1)_{B-L}$, where $h$ denotes the lightest Higgs boson with 125 GeV. Considering the constraints from the updated prediction data, we find that the production cross section of this process in the model depends on some parameters strongly.


I. INTRODUCTION
Among the particles predicted by the Standard Model (SM), the Higgs boson was the last particle discovered. Its discovery proves that correctness of the particles mass produced by the spontaneous symmetry breaking. Nevertheless, some other questions also arise. For example, is the discovered particle really the Higgs boson predicted by the SM? Are there other neutral or charged scalar fields? Whether or not the coupling of Higgs with matter fields and gauge particles meets the theoretical predictions of the SM? In addition, the SM can not provide the candidates for dark matter, can not explain the asymmetry of matter and antimatter in the universe, etc., so these require us to look for new physics beyond the SM. The B − L Symmetric SM (B-LSSM) is one of the simplest extension models of the Minimal Supersymmetric Standard Model (MSSM), which is based on the gauge symmetry group SUð3Þ C ⊗ SUð2Þ L ⊗ Uð1Þ Y ⊗ Uð1Þ B−L , where B stands for the baryon number and L stands for the lepton number, respectively.
The B-LSSM alleviates the hierarchy problem arisen in the MSSM, because the exotic singlet Higgs boson and right-handed (s)neutrinos [1][2][3][4][5][6][7] can alleviate the constraints from the experimental data of LHC, Tevatron and LEP. Furthermore, Uð1Þ B−L gauge group can help to understand the possible broken ways of R parity in the supersymmetric models [8][9][10].
In this work, we analyze the cross section of e þ e − → hhZ in the B-LSSM. The cross section and angular distribution of this process can be used to determine the self-coupling of the Higgs at future collider experiments [23]. The process has been studied in the frameworks of the SM, the two-Higgs-doublet model [13,[23][24][25], etc. In Ref. [13], the cross section in the SM was calculated as 0.16 fb at E cm ¼ 500 GeV, 0.12 fb for E cm ¼ 1000 GeV. If future experimental observations are much larger than the theoretical predictions of the SM, which can be considered as an evidence of new physics beyond the SM. Regarding the cross section, there are no experimental observations, but the production cross section for the e þ e − → hhZ process is typically of the order of 0.1 fb at the collision energy just above the threshold at about 400 GeV, and at the international linear collider with a center-of-mass energy of 500 GeV, the trilinear Higgs boson coupling can be measured via this process [26]. In addition, if the experimental measurement value deviates obviously from the SM value in the future, the model can be considered as an explanation to account for the deviations. If the experimental measured values are consistent with the SM, it will constrain our parameter space.
The paper is organized as follows. In Sec. II, we introduce the B-LSSM in detail. In Sec. III, we analyze the dependence of the cross section and angular distribution of the final state particles on the parameters in this model. The numerical results are given in Sec. IV, and some conclusions are summarized in Sec. V.

II. INTRODUCTION OF THE MODEL
In this section, we briefly introduce the basic characteristics of the B-LSSM. About the B-LSSM, there are several different versions. Here, we apply the version described in Refs. [27][28][29][30]. The B-LSSM of this version is encoded in SARAH [31][32][33][34][35], which can produce the mass matrices and interaction vertexes of this model. In this model, the local gauge group is enlarged to SUð3Þ C ⊗ SUð2Þ L ⊗ Uð1Þ Y ⊗ Uð1Þ B−L , where the Uð1Þ B−L is an additional gauge symmetry.
With the local gauge group SUð3Þ C ⊗ SUð2Þ L ⊗ Uð1Þ Y ⊗ Uð1Þ B−L , the superpotential of the B-LSSM is written as Where W MSSM denotes the superpotential of the MSSM [36], W ðB−LÞ denotes the additional terms in this model, and can be written as The quantum numbers of the superfields of the quarks and leptons are assigned aŝ The quantum numbers of chiral singlet superfields arê η 1 ∼ ð1; 1; 0; −1Þ,η 2 ∼ ð1; 1; 0; 1Þ, and that of three generations of right-handed neutrinos isν c i ∼ ð1; 1; 0; 1=2Þ. Correspondingly, the soft breaking terms in the B-LSSM are written as where the L MSSM soft denotes the soft breaking terms in the MSSM [36], and where λ B , λ 0 B represent the gauginos of Uð1Þ Y and Uð1Þ B−L , respectively. The local gauge symmetry SUð3Þ C ⊗SUð2Þ L ⊗ Uð1Þ Y ⊗Uð1Þ B−L is broken down to the electromagnetic symmetry Uð1Þ em when the Higgs fields acquires the nonzero vacuum expectation values (VEVs): Similar to the ratio of nonzero VEVs of H 1 and H 2 , we take tan β 0 ¼ u 2 u 1 denoting the ratio of nonzero VEVs of two chiral singlet superfieldsη 1 andη 2 here.
There is the gauge kinetic mixing −κ Y;BL A 0Y μ A 0μ;BL from two local Uð1Þ gauge groups, and the mixing term satisfies the gauge invariance, where A 0Y μ , A 0μ;BL represent the gauge fields of two gauge groups Uð1Þ Y and Uð1Þ B−L , respectively, and the antisymmetric tensor −κ Y;BL represents the mixing between two Uð1Þ gauge fields. The choice κ Y;BL ¼ 0 is unnatural because the mixing at low energy scale still can acquire a nonzero value through the evolution of renormalization group equations (RGEs) [37][38][39][40][41][42][43], even if we choose κ Y;BL ¼ 0 at the great uniform theory scale. The soft breaking parameters TðT ij ν ; T ij x ; T d;33 ; T u;33 Þ are proportional to the corresponding Yukawa couplings, i.e., the trilinear scalar terms in the soft supersymmetry breaking potential).
Because of the reasons above, the covariant derivative is usually written as orthogonal matrix. We can choose a proper R and then rewrite the coupling matrix as where g 1 is the hypercharge coupling constant of the SM, which can be modified in the B-LSSM. Meanwhile, two U(1) gauge fields are redefined as The gauge kinetic mixing induces some interesting phenomenology. First, A BL μ boson can mix with the A Y μ and V 3 μ bosons at the tree level. In the interaction basis (A Y μ , V 3 μ , A BL μ ), the mass squared matrix of neutral gauge bosons is written as This mass squared matrix can be diagonalized by a unitary matrix, and the mass eigenstates can be written as linear Here, θ W , θ 0 W represent two mixing angles [44]: ð13Þ ð14Þ . When x ≪ 1, the eigenvalues of Eq. (11) can be written as [45] The effective potential can be written as [46]: Here, V 0 denotes the scalar potential at tree level, ΔV 1;t represents the correction from top quark, and ΔV 2;t represents the corrections from scalar top quarks, ΔV 3;ν denotes the corrections from neutrinos, and ΔV 4;ν denotes the corrections from sneutrinos ΔV 5;b represents the correction from bottom quark, and ΔV 6;b represents the corrections from scalar bottom quarks, respectively. The concrete expressions of those pieces are masses of right-handed neutrinos, and mν iR represents the masses of right-handed sneutrinos.
The stability conditions are the detailed expression about Eq. (19) is given in Appendix A. Furthermore, the gauge kinetic mixing induces the mixing among the H 1 1 , H 2 2 ,η 1 ,η 2 at tree level. In the interaction basis (ReH 1 1 , ReH 2 2 , Reη 1 , Reη 2 ), the mass squared matrix for CP-even Higgs bosons is written as where the abbreviations are Here Δm 2 h represents the one-loop correction to mass matrix squared; the detailed expression about this is given in Appendix A, and it can be obtained by the second derivative of the effective potential. In addition, g 2 ¼ g 2 1 þ g 2 2 þ g 2 YB . The mass matrix M 2 h can be diagonalized by the 4 × 4 unitary matrix Z H .
In the interaction basis (ImH 1 1 ; ImH 2 2 ; Imη 1 ; Imη 2 ), the mass squared matrix of CP-odd Higgs can be written as: here, Δm 2 A 0 represents the one-loop correction to mass squared matrix of CP-odd Higgs, and Δm 2 A 0 can be diagonalized by the 4 × 4 unitary matrix Z A . The eigenvalues of Eq. (22) can be written as where m ð0Þ A 1;2 is the contribution under tree level approximation, and the detailed expression about θðΔmÞ is given in Appendix A. The eigenstates corresponding to Goldstone are The Higgs self-coupling can be defined as A t is the trilinear couplings between Higgs and scalar top quarks, and μ denotes the mass parameter of Higgsino. In addition, the detailed expression about tree level correction λ ð0Þ h i h j h k is given in Appendix B. The issues we discuss also involve the coupling of two CP-odd Higgs and one CP-even Higgs. The corresponding expression can be found in Appendix C.

III. CROSS SECTION OF THE HIGGS BOSON PAIR
PRODUCTION THROUGH e + e − → hhZ In this section, we will introduce the production of the Higgs boson pair through e þ e − → hhZ. The channel of the production of the Higgs boson pair is open when the collision energy E cm of the initial state particle is more than about 340 GeV. The decay and production of Higgs boson have been discussed extensively [47][48][49]. In the framework of the B-LSSM, we mainly discuss the Higgs boson pair production through e þ e − → hhZ here. We will carefully analyze the influence of relevant parameters on the total reaction cross section and angular distribution of the differential cross section in this model.
The Feynman diagrams contributing to this process are given in Figs. 1 and 2, where N denotes Z and Z 0 bosons, A represents CP-odd Higgs fields. The diagrams in Fig. 1 originate from the SM sector and the new physics sector, while those of Fig. 2 originate from the new physics sector, respectively.
In our calculation, we choose collision energy E cm ¼ 500 GeV, so we ignore the masses of positron and electron. In addition, we neglect the Feynman diagrams generated by Yukawa couplings of electron.
In the B-LSSM, additional Feynman diagrams that contribute to this process are already given in Fig with p L;R ¼ 1∓γ 5 2 , where the p 1 , p 2 , q 2 represent the momenta of the initial state particles and final state Higgs boson, k denotes the momentum of the final state vector boson, respectively. The corresponding effective amplitude can be written as: The differential cross section of this process can be written as Here, M denotes the amplitude of all of these diagrams drawn in Figs. 1 and 2, and it is written as: the Wilson coefficients of those operators are given in Appendix D. In addition, φ 1 denotes the angle between the projection of the momentum direction of the final state Higgs on the x-y plane and the x axis, Ω is the spatial solid angle between the initial state electron and the final state Higgs, and dΩ ¼ sin θ 3 dθ 3 dφ 3 , respectively. We take the momentum direction of the final state Higgs as z axis, and the momentum direction of the initial state electron on the x-z plane. Here θ 3 stands for the angle between initial state electron and final state Higgs, φ 3 is the angle between the projection of the momentum direction of the initial state electron on the x-y plane and the x axis, as shown in Fig. 3. Furthermore, E 1 , E 2 both are the energy of the final state particles Higgs, where In addition, the cross section about this process is

IV. NUMERICAL RESULTS
In this section, we will present the numerical results of the process e þ e − → hhZ. The input parameters related to the SM are chosen as m W ¼80.385GeV, m Z ¼90.1876 GeV, which constrains the parameter space of our model concretely [51]. We choose these parameters so that the corresponding theoretical prediction of the mass of the lightest CP-even Higgs fits the experimental data with 3 standard deviations: 124.37 GeV ≤ m h ≤ 125.81 GeV.
The updated experimental data [52] on searching Z 0 indicates M Z 0 ≥ 4.05 TeV at 95% Confidence Level (CL), and we choose M Z 0 ¼ 4.2 TeV in our following numerical analysis. In addition, Refs. [53,54] give us a lower bound on the ratio between the Z 0 mass and its gauge coupling at 99% CL as then the scope of g B is limited to 0 < g B ≤ 0.7. The LHC experimental data also constrains the parameter tan β 0 as tan β 0 < 1.5 [29]. In order to coincide with the constraints from the direct searches of the squarks at the LHC [55,56] and the observed Higgs signal in Ref. [57], for those parameters in the soft breaking terms, we take ; it can be obtained from Eq. (15). Now, we present our numerical results.
We adopt the latest predicted data of the SM to proceed in our analysis; it is about 0.12-0.20 fb [58]. Taking SM as a low energy effective theory, after integrating the heavy freedom of the high energy scalar, the new physical effect is constrained in some operators with dimensions greater than or equal to 6 composed of SM fields, where the operators with dimensions equal to 6 are given in Eq. (1) in the Ref. [58]. At present, the corresponding high-dimensional operator Wilson coefficients constrained by the electroweak precision test of the CMS and ATLAS experimental groups are as follows [58]: Combining the tree diagrams and the correction of these high-dimensional operators, the SM theoretical prediction on σðe þ e − → hhZÞ ≃ 0.160868 AE 0.08 fb can be obtained. Taking E cm ¼ 500 GeV, tan β 0 ¼ 1.1, mν ;33 ¼ 600 GeV, mL ;33 ¼ 600 GeV, Y ν;33 ¼ 5 × 10 −4 , and Y x;33 ¼ 0.8, we plot the total cross section σ of e þ e − → hhZ versus the parameter tan β in the Fig. 4, where the gray band represents the SM prediction with 3 standard deviations [58]. In the Fig. 4(a), we take g B ¼ 0.7, and g YB ¼ −0.9 (solid line), g YB ¼ −0.7 (dashed line), g YB ¼ −0.5 (dotted line), respectively. In the Fig. 4(b), we take g YB ¼ −0.9, and g B ¼ Obviously the theoretical prediction on the total cross section σ depends on the parameter tan β strongly. Along with the increasing of tan β, the cross section decreases steeply as tan β ≤ 15. As tan β > 20, the dependence of total cross section on tan β is mild. Furthermore, the difference between the prediction from the SM and that from the B-LSSM exceeds 3 standard deviations.
In order to further analyze how the new parameters g YB and g B in the B-LSSM affect the cross section σ, we plot the Fig. 6, and the gray band also represents the SM prediction with 3 standard deviations [58]. Taking E cm ¼ 500 GeV, tan β 0 ¼ 1.1, mν ;33 ¼ 600 GeV, mL ;33 ¼ 600 GeV, Y ν;33 ¼ 5 × 10 −4 , and Y x;33 ¼ 0.8, we plot the total cross section of e þ e − → hhZ versus the new parameters g YB and g B in the Fig. 6. In the Fig. 6(a), we take tan β ¼ 38, and g B ¼ 0.7 (solid line), g B ¼ 0.5 (dashed line), g B ¼ 0.3 (dotted line), respectively. In the Fig. 6(b), we take g YB ¼ −0.7 and tan β ¼ 35 (solid line), tan β ¼ 25 (dashed line), tan β ¼ 15 (dotted line), respectively. The total cross section σ in the B-LSSM can exceed that in the SM easily when g B and jg YB j is small. In addition, the theoretical prediction on the total cross section σ depends on the new parameters g B and g YB strongly. In the Fig. 6(a), with the decreasing of jg YB j, the total cross section increases sharply. In the Fig. 6(b), the total cross section σ decreases steeply when g B increases.
Taking E cm ¼ 500 GeV, tan β 0 ¼ 1.1, mν ;33 ¼ 600 GeV, Fig. 7, where the gray band represents the SM prediction with 3 standard deviations [58]. In the Fig. 7, we take g YB ¼ −0.9 (solid line), g YB ¼ −0.7 (dashed line), g YB ¼ −0.5 (dotted line), respectively. Obviously, the dependence of total cross section on Yukawa coupling Y x;33 is mild, and, with the decreasing of jg YB j, the total cross section increases. The mass of sneutrino in the oneloop effective potential is much smaller than the mass of the stop quark and can be almost ignored, so Y x;33 has a small effect on the total cross section σ. In addition, Y ν;33 is in the order of 10 −4 , so it also has little effect on the result.
We also plot the figure with the total cross section σ of e þ e − → hhZ versus the parameter tan β in the MSSM. In order to compare with the MSSM in Ref. [59], we adopt μ ¼ −1 TeV, E cm ¼ 500 GeV, mL ;33 ¼ 100 GeV,    8. Along with the increasing of tan β, the total cross section increases as tan β ≤ 5. When tan β > 10, the dependence of total cross section on tan β is mild. In addition, the total cross section is always within the prediction range of SM.
Finally we present the relative correction on the production cross section δ ¼ σ−σ 0 σ 0 versus the parameter tan β, where σ 0 denotes the theoretical prediction of the total cross section in the SM. Taking E cm ¼ 500 GeV, tan β 0 ¼ 1.1, g B ¼ 0.4, g YB ¼−0.4, mν ;33 ¼ 600 GeV, mL ;33 ¼ 600 GeV, Y ν;33 ¼ 5 × 10 −4 , and Y x;33 ¼ 0.8, we plot the relative correction versus tan β in the Fig. 12. The theoretical prediction on the production cross section of the B-LSSM deviates that of the SM obviously as tan β < 10. The dependence of the relative correction on tan β changes mildly as tan β > 25.

V. CONCLUSION
In this work we analyze the production cross section of e þ e − → hhZ and the self-coupling of Higgs in the B-LSSM. Some parameters affect the theoretical prediction on the production cross section of e þ e − → hhZ strongly, for example the new gauge coupling g YB . Actually the theoretical prediction on the cross section deviates from that of the SM obviously under some assumptions on the parameters of the model. Nevertheless, the correction from one-loop effective potential to the self-coupling of the lightest Higgs can be neglected safely.
Although the cross section has not been measured when the center-of-mass energy is 500 GeV, the trilinear Higgs boson coupling can be measured at international linear collider through this process. In the future, if the experimental value greatly exceeds the SM, the model can be used to explain the deviation. If the experimental value is consistent with the SM, it will constrain our parameter space. The detailed expression about stability condition: The detailed expression about Δm Here, fðQ 2 ;m 2 t;t 1;2 Þ¼ 1 4 m 2 t;t 1;2 ðln Here, the detailed expression about tree level correction λ A i A j h k is written as:   (27), respectively. Furthermore, λ h i h j h k has been given in Eq. (26) and Appendix B, and λ A i A j h k denotes the triple Higgs selfcoupling of CP-even and CP-odd Higgs; the detailed expression about this has been given in Appendix C. The superscripts (a, b, c, a 0 ; b 0 ; c 0 ; d 0 ) respectively represent the corresponding Feynman diagram labels in Figs. 1 and 2, and m h , m A denote the masses for Higgs and pseudoscalar Higgs, with i, k ¼ 1, 2, 3, 4 denoting the index of generation. ð2g YB þ g B Þ sin θ 0 W ; g eeZ 0 1 ¼ 1 2 ððg 1 sin θ W − g 2 cos θ W Þ sin θ 0 W þ ðg YB þ g B Þ cos θ 0 W Þ; g eeZ 0 2 ¼ g 1 sin θ 0 W sin θ W − 1 2 ð2g YB þ g B Þ cos θ 0 W : ðD8Þ