\subsection{Running with full top quark mass dependence ({\tt mtdep=3})}
\label{sec:appendix1}


In this appendix we give some further details about the running mode with full top quark mass dependence.

The two-loop virtual amplitudes in the NLO calculation with full top quark mass dependence are computed via a grid
which encodes the dependence of the virtual two-loop amplitude on the
kinematic invariants $\hat{s}$ and
$\hat{t}$~\cite{Heinrich:2017kxx}. We emphasize that the numerical
values $m_H=125$\,GeV and $m_t=173$\,GeV are hardcoded in this
grid and therefore should not be changed in {\tt powheg.input-save} when running in the {\tt mtdep=3} mode.
The grid is generated using python code and is directly interfaced to the
\POWHEGBOX{} fortran code via a python/C API. In order for the grid to be found by the code, the
files ({\tt events.cdf, createdgrid.py, Virt\_full\_*E*.grid}) from the
folder {\tt Virtual} need to be copied into the local folder where the
code is run. Instead of copying the files, we suggest to create a
symbolic link to the needed files.
All this is done automatically if you use the script  {\tt run.sh}.

To do this manually:
assuming the code is run from a
subfolder (e.g. {\tt testrun}) of the process folder, the link can be created in this subfolder as follows:
\begin{description}
\item{\tt ln -s ../Virtual/events.cdf events.cdf}
\item{\tt ln -s ../Virtual/creategrid.py creategrid.py}
\item{\tt for grid in ../Virtual/Virt\_full\_*E*.grid; do ln -s \$grid; done}
\end{description}
Once the links are in place, the code can be run with {\tt mtdep=3} as usual.  
The python code {\tt creategrid.py} will then combine the virtual
grids generated with the 23 combinations of coupling values
to produce a new file {\tt Virt\_full\_*E*.grid} corresponding to the
values of $c_{hhh}, c_t, c_{tt}, c_{ggh}, c_{gghh}$ defined by the user in the {\tt powheg.input-save} file.

The python code for the grid relies on the {\tt numpy} and {\tt sympy} packages, which 
the user should install separately. When building the {\tt ggHH} process the Makefile will find the embedded python 3 library
via a call to {\tt python3-config}, which the user should ensure is configured correctly and points to the correct library.
Note that on some systems the python/C API does not search for packages (such as {\tt numpy} and {\tt sympy}) in the same 
paths as the python executable would, the user should ensure that these packages can be found also by an embedded python program.
To ensure that the linked files are found, we recommend to add the run subfolder to {\tt PYTHONPATH}.

\subsection{Powheg input and run scripts}

The {\tt run.sh} script in the {\tt testrun} folder allows 
the different stages of \POWHEG{} to be run easily.
By typing {\tt ./run.sh} without any argument a menu with the
4 {\tt mtdep} running modes described above is shown. 
For all {\tt mtdep} running modes, {\tt run.sh} will make the code go through 
the various steps (parallel stages) of the calculation: 
\begin{description}
 \item[{\tt parallelstage=1}:]
generation of the importance sampling grid for the Monte Carlo integration; 
 \item[{\tt parallelstage=2}:] calculation of the integral for the inclusive cross section and an upper bounding function of the integrand;
 \item[{\tt parallelstage=3}:] upper bounding factors for the generation of radiation are computed;
 \item[{\tt parallelstage=4}:] event generation, i.e. production of {\tt pwgevents-*.lhe} files.
\end{description}

\noindent Please note: if you use the script {\tt run.sh [mtdep]}, the value for {\tt mtdep} given as an argument to {\tt run.sh} will be used, even if you specified a different value for {\tt mtdep} in {\tt powheg.input-save}.

After running {\tt parallelstage=4}, the LHE files produced by \POWHEG{} can be directly showered by either \pythia{} or \herwig{}. We provide a minimal setup for producing parton-shower matched distributions in {\tt test-pythia8}, respectively {\tt test-herwig7}. Both the angular-ordered and the dipole-shower implemented in \herwig{} can be used by changing the {\tt showeralg} flag to either {\tt default} or {\tt dipole} in {\tt HerwigRun.sh}.


Further, we should point out that 
\POWHEG{} offers the possibility to use a damping factor $h=\texttt{hdamp}$ of the
form~\cite{Alioli:2008tz,Alioli:2009je}
\begin{align}
  F=\frac{h^{2}}{(p_{T}^{hh})^2+h^{2}}\,,
\end{align}
where $p_{T}^{hh}$ is the transverse momentum of the Higgs boson pair, to
limit the amount of hard radiation which is exponentiated in the
Sudakov form factor. The setting $F\equiv1$, corresponding to {\tt hdamp}$=\infty$, results in
quite hard tails for observables like
$p_{T}^{hh}$~\cite{Heinrich:2017kxx,Heinrich:2019bkc}. Changing the damping factor $F$ by
setting the flag {\tt hdamp} to some finite value in the input
card softens the high transverse momentum tails. Varying {\tt hdamp} allows shower uncertainties to be assessed
within the \POWHEG{} matching scheme. However, 
{\tt hdamp} should not be so low that it starts to cut into the Sudakov
regime. In fact, a too low value for {\tt hdamp} could spoil the
logarithmic accuracy of the prediction. For this reason we suggest not
to choose values for {\tt hdamp} below $\sim 200$. Our default value is  {\tt hdamp=250}.


\section{Introduction}

\input{intro.tex}

\section{Anomalous couplings in Higgs boson pair production}
\label{sec:couplings}

\input{couplings.tex}

\section{Description of the code}
\label{sec:usage}

\input{usage.tex}

\section{Phenomenological results}
\label{sec:results}

\input{results.tex}


\clearpage

\section{Conclusions}
\label{sec:conclusions}

We have presented a publicly available implementation of Higgs boson pair production in gluon
fusion within an Effective Field Theory framework, calculated at full NLO QCD, in the \POWHEGBOXV{}.
The code allows five anomalous couplings relevant for di-Higgs
production to be varied and offers the possibility to produce fully differential final states.

We have also investigated the behaviour of the shape of the invariant
mass distribution of the Higgs boson pair, $m_{hh}$, for seven benchmark points
representing characteristic $m_{hh}$-shapes based on an NLO analysis~\cite{Capozi:2019xsi}.
In addition, for one of the benchmark points (benchmark 4), characterised by an enhanced tail in the $m_{hh}$ distribution,
we carried out a comparison to the full NLO SM result as well as to the $m_t\to \infty$ approximation, including scale uncertainties, for both the $m_{hh}$ and the $p_{T}^{hh}$ distributions.
We found that the two distributions show different characteristics concerning the distinction of the BSM curve from the SM (and approximate SM) curves: 
while in the  $m_{hh}$ distribution the enhanced tail of benchmark 4 is clearly outside the uncertainty bands of the SM predictions, in the $p_{T}^{hh}$ distribution the bands for benchmark 4 and the SM in the $m_t\to \infty$ approximation overlap. This again demonstrates the importance of using full NLO predictions, particularly in trying to resolve partly degenerate directions in the space of anomalous couplings.


We further produced results for the  $p_{T}^{hh}$ distribution  matched to three different parton showers: \pythia{} and two different \herwig{} parton showers. The two \herwig{} parton showers show very similar results.
From previous SM results, it is known that loop-induced processes like $gg\to HH$ in \POWHEG{} matched to \pythia{} can show a harder tail than with other parton showers.
However, in the case of benchmark point 4,  \pythia{} produces less additional hard radiation than in the SM case, such that the \pythia{}  and \herwig{} results are much more similar.

Our studies show that the behaviour of higher-order effects known from the SM does not necessarily carry over to the BSM
case, such that precise predictions for both cases are necessary to clearly identify new physics effects.
With our code, fully exclusive studies of anomalous couplings in Higgs boson pair production at NLO QCD are possible.


\section*{Acknowledgements}
We would like to thank Gerhard Buchalla for useful discussions and
Matteo Capozi and Gionata Luisoni for collaboration on earlier
versions of this code. We also are grateful to Tom\'a\v{s} Je\v{z}o
and Emanuele Re for very useful comments about the \POWHEGBOX.
This research was supported in part by the COST Action CA16201 (`Particleface') of the European Union
and  by  the  Deutsche  Forschungsgemeinschaft (DFG, German Research
Foundation) under grant 396021762 - TRR 257. LS is supported by the Royal Society under 
grant number RP\textbackslash R1\textbackslash 180112.
MK acknowledges support by the Swiss National  Science  Foundation  (SNF)  under  grant  number  200020-175595
and by the Forschungskredit of the University of Zurich, grant no. FK-19-102.
We also acknow\-ledge
resources provided by the Max Planck Computing and Data
Facility (MPCDF).


\renewcommand \thesection{\Alph{section}}
\renewcommand{\theequation}{\Alph{section}.\arabic{equation}}
\setcounter{section}{0}
\setcounter{equation}{0}

\section{Appendix}

\input{appendix.tex}


\bibliographystyle{JHEP}




\subsection{Structure of the code}

The code is an extension of the one presented in Ref.~\cite{Heinrich:2019bkc} to
include the possibility of varying all five anomalous couplings rather
than only the trilinear Higgs coupling.
For the virtual two-loop corrections,  we have built on the results of the
calculations presented in Refs.~\cite{Borowka:2016ehy,Borowka:2016ypz}.
These results were obtained by performing a (partial) reduction of the two-loop amplitude to master integrals based on {\sc Reduze}~\cite{Studerus:2009ye,vonManteuffel:2012np} and a subsequent numerical evaluation of the master integrals using the program {\sc SecDec}~\cite{Borowka:2015mxa,Borowka:2017idc}.
This has been done with the top quark- and Higgs masses fixed to a numerical value. Therefore these mass values should not be changed in the {\tt ggHH} code.

The real radiation matrix elements were implemented
using the interface between \textsc{GoSam}{}~\cite{Cullen:2011ac,Cullen:2014yla} and
the \POWHEGBOX~\cite{Alioli:2010xd,Luisoni:2013cuh}.
The extra matrix elements occurring in the EFT framework have been
generated by \textsc{GoSam}{}{} via a model file in {\tt
  UFO} format~\cite{Degrande:2011ua} which has
been developed in Ref.~\cite{Buchalla:2018yce}, 
derived from the effective Lagrangian in Eq.~(\ref{eq:ewchl}) using {\sc FeynRules}~\cite{Alloul:2013bka}.

The framework presented in Ref.~\cite{Heinrich:2019bkc} to interface the two-loop 
virtual contribution in {\tt POWHEG} is generalised in the following way: instead 
of a second-order polynomial (as for variations of $c_{hhh}$ only), at NLO we can write the 
squared matrix element for variations of all five anomalous couplings as in Eq.~(\ref{eq:AcoeffsNLO}),
following Refs.~\cite{Azatov:2015oxa,Carvalho:2015ttv,Buchalla:2018yce}.

\begin{align}
|\mathcal{M}_{\rm BSM}|^2  &=   a_1\, c_t^4 + a_2 \, c_{tt}^2  + a_3\,  c_t^2 c_{hhh}^2  +
a_4 \, c_{ggh}^2 c_{hhh}^2  + a_5\,  c_{gghh}^2  +
a_6\, c_{tt} c_t^2 + a_7\,  c_t^3 c_{hhh} \nonumber\\
& + a_8\,  c_{tt} c_t\, c_{hhh}  + a_9\, c_{tt} c_{ggh} c_{hhh} + a_{10}\, c_{tt} c_{gghh} +
a_{11}\,  c_t^2 c_{ggh} c_{hhh} + a_{12}\, c_t^2 c_{gghh} \nonumber\\
& + a_{13}\, c_t c_{hhh}^2 c_{ggh}  + a_{14}\, c_t c_{hhh} c_{gghh} +
a_{15}\, c_{ggh} c_{hhh} c_{gghh}\, + a_{16}\, c^3_t c_{ggh} \nonumber\\
& + a_{17}\,  c_t c_{tt} c_{ggh}
+ a_{18}\, c_t c_{ggh}^2 c_{hhh} + a_{19}\, c_t c_{ggh} c_{gghh} \, + a_{20}\,  c_t^2 c_{ggh}^2
\nonumber\\
&+ a_{21}\, c_{tt} c_{ggh}^2
+ a_{22}\, c_{ggh}^3 c_{hhh} + a_{23}\, c_{ggh}^2 c_{gghh}  \,.
\label{eq:AcoeffsNLO}
\end{align}

For the Born-virtual interference term, we produce grids using 6715 points (5194 points at $\sqrt{s}=14$ TeV and 1521 points at $100$ TeV)
for 23 linearly independent sets of couplings. This enables us to 
derive, for each phase-space point, the coefficients $a_1, \dots, a_{23}$ by 
interpolation. Once the user has chosen a set of anomalous couplings, the 23 grids are
combined into one using Eq.~(\ref{eq:AcoeffsNLO}). This step is performed only once, in 
the first {\tt POWHEG} parallel stage. The Born and real contributions are evaluated exactly
for the chosen anomalous couplings without relying on a grid or interpolation. Note that
the $a_i$ coefficients of Eq.~(\ref{eq:AcoeffsNLO}) are not equal to the $A_i$ coefficients of 
Ref.~\cite{Buchalla:2018yce}, which are derived for the (normalised) cross-section and 
not for the Born-virtual interference term.

The original phase space points were produced with SM couplings,
therefore some regions, for example the low $m_{hh}$-region, are less
populated than the peak of the $m_{hh}$-distribution in the SM.
The statistical uncertainty on the input data induces a systematic uncertainty
on the 23 grids.
We have checked the relative size of the uncertainties of the virtual
corrections in each bin of the $m_{hh}$ distributions for all seven
benchmark points. This uncertainty is below 2\% throughout the
whole $m_{hh}$ range, except for the first bin.
This bin is poorly populated in the SM and therefore the uncertainties
in this bin are larger for coupling configurations where the low-$m_{hh}$
region is very different from the SM case or where the relative size of the virtual contribution is large. 
Hence we find uncertainties in the first $m_{hh}$ bin of about 6\%  for all benchmarks except 
for benchmark 5, which has a 12\% uncertainty in this bin.


\subsection{Usage of the code}

The code  can be found at the web page
\begin{center}
 \url{http://powhegbox.mib.infn.it}
\end{center}
under {\tt User-Processes-V2} in the {\tt ggHH} process directory.
An example input card ({\tt powheg.input-save}) and a run script ({\tt run.sh}) 
are provided in the {\tt testrun} folder accompanying the code.

In the following we only describe the input parameters that are specific to the
process $gg\to HH$ including five anomalous couplings. The
parameters that are common to all \POWHEGBOX{} processes can be found
in the \POWHEG{} manual {\tt V2-paper.pdf} in the \POWHEGBOXV{}/{\tt Docs} directory.


\subsubsection*{Running modes}
The code contains the SM NLO QCD amplitudes with full top quark mass dependence.
A detailed description of the different approximations can be found in
Ref.~\cite{Borowka:2016ypz}. For the Standard Model case as well as for BSM-values of the trilinear
Higgs coupling $c_{hhh}$, the code can be run in four different modes, either by
changing the flag {\tt mtdep} in the \POWHEGBOX{} run card {\tt
  powheg.input-save}, or by using the script {\tt run.sh [mtdep
  mode]}.
If all five anomalous couplings are varied, there is only the
possibility of either calculating at

\begin{itemize}
\item NLO with full top quark mass dependence, or
\item LO (setting 
{\tt bornonly=1}) in either the full theory or in the $m_t\to\infty$ limit.
\end{itemize}

In more detail, the following choices are available:
\begin{description}
 \item[{\tt mtdep=0}:]{all amplitudes
   are computed in the $m_t\to\infty$ limit (HTL).} {\em This option is only available at NLO in the SM
   case or if only $c_{hhh}$ is varied, or at LO.}
 \item[{\tt mtdep=1}:]{computation using Born-improved HTL. In this
   approximation the fixed-order part is computed at NLO in the heavy top limit
   and reweighted pointwise in the phase-space by the LO matrix
   element with full mass dependence divided by the LO matrix
   element in the HTL.} {\em This option is only available at NLO in the SM
   case or if only $c_{hhh}$ is varied, or at LO.}
 \item[{\tt mtdep=2}:]{computation in the approximation FT$_{\mathrm{approx}}$. In
   this approximation the matrix elements for the Born and the real
   radiation contributions are computed with full top quark mass dependence, whereas the virtual part is
   computed as in the Born-improved HTL. }{\em This option is only available at NLO in the SM
   case or if only $c_{hhh}$ is varied, or at LO.}
 \item[{\tt mtdep=3}:]{NLO computation with full top quark mass dependence.}
\end{description}



\subsubsection*{Input parameters}

The bottom quark is considered massless in all four {\tt mtdep} modes. The Higgs
bosons are generated on-shell with zero width. A decay can be attached
through the parton shower in the narrow-width approximation. However,
the decay is by default switched off (see the {\tt hdecaymode} flag in the
example {\tt powheg.input-save} input card in {\tt testrun}).

The masses of the Higgs boson and the top quark are set by default to
$m_h=125$\,GeV, $m_t=173$\,GeV, respectively, and their widths
have been set to zero. The full SM two-loop virtual contribution has
been computed with these mass values hardcoded, therefore they should not be changed when running with {\tt mtdep = 3}, otherwise the two-loop virtual part would contain a different top or Higgs mass from the rest of the calculation.
 It is no problem to change the values of $m_h$
and $m_t$ via the {\tt powheg.input-save} input card when running with
{\tt mtdep} set to $0$, $1$ or $2$.

The Higgs couplings as defined in the context of the Electroweak Chiral Lagrangian (see~\cite{Buchalla:2018yce} and references within), can be varied directly in the {\tt powheg.input} card. These are, with their SM values as default:
\begin{description}[leftmargin=!,labelwidth=\widthof{{\tt cgghh=0.0}:}]
\item[{\tt chhh=1.0}:] { the ratio of the Higgs trilinear coupling to its SM value,}
\item[{\tt ct=1.0}:] { the ratio of the Higgs Yukawa coupling to the top quark to its SM value,}
\item[{\tt ctt=0.0}:] { the effective coupling of two Higgs bosons to a top quark pair,}
\item[{\tt cggh=0.0}:] { the effective coupling of two gluons to the Higgs boson,}
\item[{\tt cgghh=0.0}:] { the effective coupling of two gluons to two Higgs bosons.}
\end{description}

 The possibility of varying all Higgs couplings (rather than $c_{hhh}$ only)  is only
 available in the mode {\tt mtdep=3} (full NLO).
 More details about the  {\tt mtdep=3} running mode are given in
 Appendix \ref{sec:appendix1}.

 The runtimes are dominated by the evaluation of the real radiation
 part. When run in the full NLO mode, the runtimes we observed were in
 the ballpark of 100 CPU hrs for an uncertainty of about 0.1\% on the
 total cross section.

