The power counting also changes as now we have $|\frs|=5$ and \begin{tikzpicture}[baseline={(0,-0.1)}]
\draw[kerneldiff] (0,0) -- (1,0);
\draw[rho] (-0.5,0) -- (0,0) node[int]{};
\end{tikzpicture}
now admits the bound
\begin{tikzpicture}[baseline={(0,-0.1)}]
\draw (-0.5,0) -- (1,0);
\draw (0.25,0) node[labl] {\tiny $-3,0$};
\end{tikzpicture}.
Consequently, we can bound both terms in \eqref{eq:model building 3} using Lemma \ref{lem:kernel-convolution} (ii) and (iv) in the first inequality and (i) in the second one:
\begin{equ}
\begin{tikzpicture}[scale=0.9,baseline={(0,0.8)}]
\draw (-1,1) -- (0,0) -- (1,1) node[int] {} -- (0,2) -- (-1,1) node[int] {};
\draw (0,0) node[vari]{} -- (0,2) node[int]{} ;
\draw (0,1) node[labl] {\tiny $0,-3$};
\draw (-0.6,0.5) node[labl] {\tiny $-3,0$};
\draw (-0.6,1.5) node[labl] {\tiny $0,-3$};
\draw (0.6,1.5) node[labl] {\tiny $0,-3$};
\draw (0.6,0.5) node[labl] {\tiny $0,-3$};
\end{tikzpicture}
\quad\lesssim
\begin{tikzpicture}[scale=0.9,baseline={(0,0.8)}]
\draw (0,0) to[out=45,in=-45,distance=1.5cm] (0,2);
\draw (0,0) to[out=135,in=-135,distance=1.5cm] (0,2) node[int]{};
\draw (0,0) node[vari]{} -- (0,2) node[int]{} ;
\draw (0,1) node[labl] {\tiny $0,-3$};
\draw (-0.6,0.5) node[labl] {\tiny $-1,0$};
\draw (0.6,0.5) node[labl] {\tiny $0,-1$};
\end{tikzpicture}
\quad\lesssim
\begin{tikzpicture}[scale=0.9,baseline={(0,0.8)}]
\draw (0,0) node[vari]{} -- (0,2) node[int]{} ;
\draw (0,1) node[labl] {\tiny $-1,-4$};
\end{tikzpicture}
\quad\lesssim
\big|{\log|x|_\d}\big|.
\end{equ}
Concerning the non-recentered term in the first chaos component of $\<Phi-5new>$, it can be written as
\begin{equ}
\cW^{(\eps;1)}_2\<Phi-5new>=6\,
\begin{tikzpicture}[scale=0.7,baseline={(0,0.3)}]
\draw[kerneldiff] (-1,1) -- (0,0);
\draw[kernels3] (0,0) -- (1,1) ;
\draw (-1,2.5) -- (0,1.5) -- (1,2.5) ;
\draw (0,1.5) -- (0,3) node [noise]{};
\draw[rho] (-1,1) node[int]{} -- (-1,1.75);
\draw[rho] (1,1) node[int]{} -- (1,1.75);
\draw[rho] (-1,2.5) node[int]{} -- (-1,1.75);
\draw[rho] (1,2.5) node[int]{} -- (1,1.75);
\draw (0,0) -- (0,1.5);
\draw (0,0) node[vari]{};
\draw (0,1.5) node[int]{};
\draw (-1,1.75) node[int]{};
\draw (1,1.75) node[int]{};
\end{tikzpicture}
\quad + \quad
\draw[kerneldiff] (-1,1) -- (0,0);
\draw (0,0) -- (1,1) ;
\draw (-1,2.5) -- (0,1.5) -- (1,2.5) ;
\draw (0,1.5) -- (0,3) node [noise]{};
\draw[rho] (-1,1) node[int]{} -- (-1,1.75);
\draw[rho] (1,1) node[int]{} -- (1,1.75);
\draw[rho] (-1,2.5) node[int]{} -- (-1,1.75);
\draw[rho] (1,2.5) node[int]{} -- (1,1.75);
\draw (0,0) -- (0,1.5);
\draw (0,0) node[vari]{};
\draw (0,1.5) node[int]{};
\draw (-1,1.75) node[int]{};
\draw (1,1.75) node[int]{};
\end{tikzpicture}
\quad.
\end{equ}
Therefore we have (using two instances of
\begin{tikzpicture}[baseline={(0,-0.1)}]
\draw (0,0) node[vari]{};
\end{tikzpicture}
to denote the two different running variables $x$ and $\bar x$)
by applying Lemma \ref{lem:kernel-convolution} (ii) and (iv) in the second inequality, (i) in the third one, and (iii) in the fourth and fifth one
\begin{equs}\label{eq:graph-calc-2}
\big|\scal{\cW^{(\eps;1)}_0\<Phi-5new>\,(x;\cdot),\cW^{(\eps;1)}_0\<Phi-5new>\,(\bar x;\cdot)}\big|  &  \lesssim\,
\begin{tikzpicture}[scale=0.9,baseline={(0,0.8)}]
\draw (-1,1) -- (0,0) -- (1,1) node[int] {} -- (0,2) -- (-1,1) node[int] {};
\draw (0,0) node[vari]{} -- (0,2) node[int]{} ;
\draw (0,1) node[labl] {\tiny $0,-3$};
\draw (-0.6,0.5) node[labl] {\tiny $-3,0$};
\draw (-0.6,1.5) node[labl] {\tiny $0,-3$};
\draw (0.6,1.5) node[labl] {\tiny $0,-3$};
\draw (0.6,0.5) node[labl] {\tiny $0,-3$};
\draw (2,1) -- (3,0) -- (4,1) node[int] {} -- (3,2) -- (2,1) node[int] {};
\draw (0,2) -- (1.5,2) node[int] {} -- (3,2);
\draw (0.75,2) node[labl] {\tiny $0,-3$};
\draw (2.25,2) node[labl] {\tiny $0,-3$};
\draw (3,0) node[vari]{} -- (3,2) node[int]{} ;
\draw (3,1) node[labl] {\tiny $0,-3$};
\draw (2.4,0.5) node[labl] {\tiny $-3,0$};
\draw (2.4,1.5) node[labl] {\tiny $0,-3$};
\draw (3.6,1.5) node[labl] {\tiny $0,-3$};
\draw (3.6,0.5) node[labl] {\tiny $0,-3$};
\end{tikzpicture}
\\&\lesssim\begin{tikzpicture}[scale=0.9,baseline={(0,0.8)}]
\draw (0,0) to[out=45,in=-45,distance=1.5cm] (0,2);
\draw (0,0) to[out=135,in=-135,distance=1.5cm] (0,2) node[int]{};
\draw (0,0) node[vari]{} -- (0,2) node[int]{} ;
\draw (0,1) node[labl] {\tiny $0,-3$};
\draw (-0.6,0.5) node[labl] {\tiny $-1,0$};
\draw (0.6,0.5) node[labl] {\tiny $0,-1$};
\draw (0,2) -- (3,2);
\draw (1.5,2) node[labl] {\tiny $0,-1$};
\draw (3,0) to[out=45,in=-45,distance=1.5cm] (3,2);
\draw (3,0) to[out=135,in=-135,distance=1.5cm] (3,2) node[int]{};
\draw (3,0) node[vari]{} -- (3,2) node[int]{} ;
\draw (3,1) node[labl] {\tiny $0,-3$};
\draw (2.4,0.5) node[labl] {\tiny $-1,0$};
\draw (3.6,0.5) node[labl] {\tiny $0,-1$};
\end{tikzpicture}
\\
\quad
\begin{tikzpicture}[scale=0.9,baseline={(0,0.4)}]
\draw (0,0) node[vari]{} -- (0,1) node[int]{} ;
\draw (0,0.5) node[labl] {\tiny $-1+\eta,-4$};
\draw (0,1) -- (3,1);
\draw (1.5,1) node[labl] {\tiny $0,-1$};
\draw (3,0) node[vari]{} -- (3,1) node[int]{} ;
\draw (3,0.5) node[labl] {\tiny $-1+\eta,-4$};
\end{tikzpicture}
\\
&\lesssim |x|_\d^{\eta}\quad
\begin{tikzpicture}[scale=0.9,baseline={(0,-0.1)}]
\draw (0,0) node[vari]{} -- (2,0) node[int]{} -- (4,0) node[vari]{};
\draw (1,0) node[labl] {\tiny $0,-1$};
\draw (3,0) node[labl] {\tiny $-1+\eta,-4$};
\end{tikzpicture}
\\
&\lesssim |x|_\d^{\eta}|\bar x|_\d^\eta\quad
\begin{tikzpicture}[scale=0.9,baseline={(0,-0.1)}]
\draw (0,0) node[vari]{} -- (2,0) node[vari]{};
\draw (1,0) node[labl] {\tiny $0,-1$};
\end{tikzpicture}
\\
&=|x|_\d^{\eta}|\bar x|^\eta_\d\|x-\bar x\|^{-1}\;.
\end{equs}
Here $\eta<0$ is arbitrary and introduced only to avoid equalities in the conditions of Lemma~\ref{lem:kernel-convolution}.
Since $|x|_\d\lesssim\|x-z\|$, $|\bar x|_\d\lesssim\|\bar x-z\|$ can be assumed, this proves \eqref{eq:model-conv3}. The remaining arguments are then exactly as in the proof of Lemma \ref{lem:PAM model convergence}.
\end{proof}

\section{Proofs of the main results}\label{sec:proofs}

The convergence results are obtained via \cite[Thm.~5.6]{GH17}
(or a slight modification of it),
the notation of which we shall use without repeating them here.
Let us however summarise the step which is the main difference
to most of the literature and the source of the boundary renormalisation.
This concerns setting up the convolution operator of a fixed point problem.
Given a modelled distribution $f$ -- think of the right-hand side of an abstract equation -- we will define an abstract convolution $\scP_\eps(f),$
where the operator $\scP_\eps$ is built from the following ingredients:
\begin{claim}
\item A kernel $\bar \cK$ that is $2$-smoothing in the sense of \cite[Ass.~5.1]{H0} -- this will always be the truncated (translation invariant) heat kernel for us;
\item A model $(\Pi^\eps,\Gamma^\eps)$ admissible for $\bar\cK$ -- we will always choose from the models discussed above;
\item A remainder kernel $Z_\eps$ that is smooth locally, but not necessarily globally up to the boundary;
\item A distribution $\zeta$ that agrees with $\cR^\eps f$ away from the spatial boundary $\R\times\d D$.
\end{claim}
Given these ingredients, \cite[Lem~4.12 \& 4.16]{GH17} yields a modelled distribution $\scP_\eps(f)$ that satisfies
\begin{equ}
\cR^\eps\big(\scP_\eps(f)\big)(x)=\scal{(K+Z_\eps)(x,\cdot),\zeta}.
\end{equ}
Moreover, $\scP_\eps$ is locally Lipschitz continuous in the natural topologies.
In particular, to get convergence results, one has to ensure that the choice of $\zeta$ is sufficiently stable in the $\eps\to 0$ limit.
This is far from obvious:
for the type of weighted spaces we wish to use, the reconstruction theorem
\cite[Thm~3.10]{H0} can only be applied locally in the interior of
$\R\times D$, which therefore yields a distribution $\cR^\eps f$ only on the interior.
Extending it to a distribution $\zeta$ \emph{on the whole space},
may not be obvious or indeed even unique.
Of course as long as $\eps>0$, there is trivial canonical extension, but it is \emph{not} true in general
that this choice will converge to a limit.
To ensure convergence, we modify the canonical choice by certain Dirac masses on the boundary, which in turn has the effect of changing the boundary conditions (see Section \ref{sec:lin} or \cite[Rem~1.5]{HP19}).
In line with [8], we will also use the notation $\zeta=\hat\cR^\eps f$, to emphasise its relation to the `usual' reconstruction operator from [10].

%At this point we were of course rather loose about the exact conditions on weights and regularities, we will be more precise in our concrete applications below.
Let us recall that in the weighted spaces of modelled distribution
used in [8]
weights are given by triples $w=(\eta,\sigma,\mu)\in\R^3$ describing singularities near the temporal and spatial boundaries and their intersection, respectively.
It will often be convenient to use the notation $(\eta)_3=(\eta,\eta,\eta)$.
The lowest homogeneity symbol appearing in the target space of a given space of modelled distributions will usually be denoted by $\alpha$.
There are two general situations when the issue of choosing $\zeta$ can be avoided.
First, as long as $\sigma\wedge\alpha>-1$, there \emph{is} a canonical choice for $\hat\cR$, see \cite[Thm~4.9]{GH17}, which furthermore agrees with the natural extension for smooth models.
Second, if $\sigma\wedge\alpha>-2$ and $\bar\cK+Z_\eps$ is a Dirichlet kernel, then $\scP_\eps (f)$ can be defined \emph{without} defining $\hat\cR f$, see Lemma \ref{lem:Dirichlet-convolution} below.

Finally, we address the above mentioned `slight modification' of \cite[Thm~5.6]{GH17} that we will use in the case of $\Phi^4_3$,
namely that we do not consider the remainder part of the kernel fixed
as $\eps$ varies.
It is straightforward to check that the integration operation against smooth remainders with singularities at the boundary \cite[Lem~4.16]{GH17} is Lipschitz continuous with respect to the natural norm $\|\cdot\|_{\scZ_{\beta,P}}$ on the space of remainder kernels.
One then gets the following simple extension of \cite[Thm.~5.6]{GH17}.
\begin{proposition}\label{prop:tiny extension}
Assume the setting of \cite[Thm.~5.6]{GH17},
with however allowing for $Z\in\scZ_{\beta,P}$ to depend on $\eps\in[0,1]$.
Suppose that
\begin{equ}
\|Z^\eps-Z^0\|_{\scZ_{\beta,P}}\rightarrow0.
\end{equ}
Then the conclusions of \cite[Thm.~5.6]{GH17} hold.
\end{proposition}

\subsection{PAM}
As a warm-up example of the above setup without boundary renormalisation, define
the modelled distribution
\begin{equ}
\bY_\eps=-\scP_\eps(\bone_D\<0>),
\end{equ}
with $Z_\eps=Z$ being the Neumann remainder (see either \cite[Ex~4.15]{GH17}, \cite[App~A]{HP19}, or Section~\ref{sec:robin} above) and $\zeta=\bone_D\xi_\eps$.
Here and below $\bone_D$ is the indicator function of $D$.
Note that one can easily make sense (probabilistically) of the product $\bone_D\xi$,
and one moreover has $\bone_D\xi_\eps\rightarrow\bone_D\xi$ in $\cC^{-3/2-\kappa}$.
It is also easy to see that $\bone_D\<0>\in\cD_P^{\kappa,(-3/2-\kappa)_3}$.
setting $\bw_\eps=\scD\bY_\eps$,
we have that $\bw_\eps \to\bw_0\in \cD_P^{1+\kappa,(-1/2-\kappa)_3}$ and $\cR^\eps\bw_\eps=\Psi_\eps$, for all $\eps\in[0,1]$.

This allows us to set up the abstract fixed point problem corresponding to \eqref{e:etPAM} as follows.
Let $\bone_D^+$ denote the indicator function of $\R_+\times D$.
Let $Z_\eps=Z$ be the either the Dirichlet remainder (for Theorem \ref{thm:PAM-D})
%or the Robin remainder with $a=a_\rho-\tfrac{|{\log \eps}|}{8\pi}$
or the Neumann remainder (for Theorem \ref{thm:PAM-N}).
%Here, the constant $a_\rho$ is obtained from Lemma \ref{lem:reno Psi^2 PAM}.
Let furthermore $\bi_\eps$ be the modelled distribution
responsible for the contribution of the initial condition:
it is simply obtained by taking the Taylor polynomial lift of
the function $t,x\mapsto\scal{(\bar\cK+Z)_t(x,\cdot),u_0e^{Y_\eps}}$.
The abstract counterpart of \eqref{e:etPAM} then reads as
\begin{equ}\label{eq:PAM-abstract}
\bv_\eps=\scP_\eps\big(\bone_D^+ \bv_\eps|\bw_\eps|^2\big)
\end{equ}
The operators $\scP_\eps$ and $\scP_\eps'$ will only differ in the choice of $\hat\cR^\eps$. In fact as we will shortly see, for the latter this choice will be trivial.

We now specify the spaces of modelled distributions on which \eqref{eq:PAM-abstract} can be shown to be a contraction and therefore to have a unique solution.
Denote by $T_0^{\geq \alpha}$ the subspace of $T_0$ generated by basis vectors of homogeneity greater or equal to $\alpha$
(note that no new symbols are included in these spaces)
and by $\bar T$ the polynomial subspace.
Fix the exponents $\gamma=3/2+3\kappa$, $\gamma_0=\kappa$,
$\alpha=0$, $\alpha_0=-1-2\kappa$, $\alpha_1=-1/2-3\kappa$, the sectors
\begin{equ}
V=\cI(T_0^{\geq \alpha_0})+\bar T,\quad V_0=T_0^{\geq \alpha_0},\quad V_1=T_0^{\geq \alpha_1},
\end{equ}
as well as
\begin{equs}
%\alpha &=0, &\quad
\eta &=-1/2+3\kappa&\quad \sigma &=1-3\kappa&\quad \mu &=-1/2+3\kappa;\\
%\alpha_0 &=-1-2\kappa, &\quad
\eta_0 &=-3/2+\kappa,&\quad \sigma_0 &=-1-2\kappa,&\quad \mu_0 &=-3/2+\kappa;\\
%\alpha_1 &=-1/2-3\kappa, &\quad
\eta_1 &=-2+\kappa,&\quad \sigma_1 &=-1/2-4\kappa,&\quad \mu_1 &=-2+2\kappa.
\end{equs}
We use the shorthand $w=(\eta,\sigma,\mu)$ and similarly for $w_1$ and $w_2$.

It is easy to verify that, by \cite[Lem~5.4]{GH17}, $\bi_\eps\to\bi_0$ in $\cD^{\gamma,w}(V)$, provided $\kappa$ is sufficiently small.
It also follows from \cite[Lem.~4.3]{GH17} that the mappings
\begin{equ}
f\mapsto f|\bw_\eps|^2,\quad\quad f\mapsto \scD f\cdot\bw_\eps
\end{equ}
are globally Lipschitz continuous from $\cD^{\gamma,w}(V)$
to $\cD^{\gamma_0,w_0}(V_0)$ and
$\cD^{\gamma_0,w_1}(V_1)$, respectively.
Furthermore, since $\sigma_1\wedge\alpha_1>-1$,
\cite[Thm~4.9]{GH17} provides a continuous extension $\hat \cR^\eps$ of the reconstruction operator on $\cD^{\gamma_0,w_1}(V_1)$,
which agrees with the canonical extension for $\eps>0$.
This now completes the definition of $\scP_\eps'$ in \eqref{eq:PAM-abstract}.

Furthermore, in the case of Dirichlet boundary conditions, $\scP_\eps$ is also well-defined, thanks to Lemma \ref{lem:Dirichlet-convolution}.
By \cite[Thm~5.6]{GH17} one can conclude that,
provided $\kappa>0$ is sufficiently small, it has a unique local in time solution $\bv_\eps$.
Furthermore, since the right-hand side is globally Lipschitz continuous (in fact, linear), one can easily extend the local solutions to the whole time interval $[0,1]$.
Now the proof can be easily concluded.

\begin{proof}[Proof of Theorem \ref{thm:PAM-D}]
Invoking the usual arguments concerning the renormalisation in the bulk (see e.g. \cite[Sec~9]{H0}),
one sees that the function $v_\eps:=\cR^\eps\bv_\eps$ solves \eqref{e:etPAM} with $0$ Dirichlet boundary conditions.
Hence $v_\eps e^{-Y_\eps}$ solves solves \eqref{e:ePAM},
trivially also with $0$ Dirichlet boundary conditions,
and therefore it coincides with $u_\eps^{\Dir}$.
By \cite[Thm.~5.6]{GH17}, $\bv_\eps\to\bv_0$, the latter
clearly does not depend on $\rho$,
and so the limit $\lim_{\eps\to0}u_\eps^{\Dir}=(\cR^0\bv_0)e^{-Y_0}=:u^{\Dir}$ is also independent of $\rho$, yielding the claim.
\end{proof}

In the Neumann case, there is no easy way to construct (or avoid constructing, as above) $\hat\cR^\eps$ for the integration operator $\scP_\eps$.
The main remaining step is therefore to build $\hat\cR^\eps$ as a continuous map from
$\cV_\eps:=\cD^{\gamma,w}(V)|\bw_\eps|^2$ to $\CC^{-3/2+\kappa}$,
which is where all previous steps in the article will be used.
First of all, notice that $\bw_\eps$ is of the form
\begin{equ}\label{eq:qwe}
\end{equ}
for some $(\bar T)^3$-valued $\varphi_\eps$.
One however can use the newly added symbols to out regularity structure to write an equivalent function
\begin{equ}\label{eq:rty}
\end{equ}
where $\tilde\varphi_\eps(x)=\E|\Psi_\eps|^2(x)-\ell_\eps(\<PAM-2a-Small>)$.
Although the latter expression seems ill-defined for $\eps=0$,
it is clear that $\tilde\varphi_\eps$ has pointwise limits away from the boundary.
In fact, in the notation of Lemma \ref{lem:reno Psi^2 PAM} we have the expression $\tilde \varphi_\eps(x)=R^2_\eps+R^3_\eps$.
The $|\cdot|^2$ in \eqref{eq:rty} is merely notational and serves only as an analogy with \eqref{eq:qwe}.
One then trivially has $\cR^\eps|\bw_\eps|^2=\cR^\eps|\tilde \bw_\eps|^2$,
in fact one has $\cR^\eps f|\bw_\eps|^2=\cR^\eps f|\tilde\bw_\eps|^2$
for all $f\in\cD^{\gamma,w}(V)$.
Let us therefore define the isometry $\iota_\eps:\cV_\eps\to\tilde \cV_\eps:=\cD^{\gamma,w}(V)|\tilde\bw_\eps|^2$ by replacing the $|\bw_\eps|^2$ factor with $|\tilde\bw_\eps|^2$.
We can conclude that it suffices to construct $\hat \cR^\eps$ on $\tilde \cV_\eps$, and then it can be pulled back to $\cV_\eps$ by $\iota_\eps$.

\begin{lemma}\label{lem:hatR-PAM}
For all $\eps\in[0,1]$ there exist maps $\hat \cR^\eps:\tilde\cV^\eps\to\CC^{-3/2+\kappa}$ with the following properties.
\begin{claim}
\item For all $\eps\in[0,1]$, $g\in\tilde\cV^\eps$, and $\psi\in\cC^\infty_0(\R_+\times D)$ one has
\begin{equ}\label{eq:hatR-PAM1}
\scal{\hat\cR^\eps g,\psi}=\scal{\cR^\eps g ,\psi};
\end{equ}
\item The bound
\begin{equ}\label{eq:hatR-PAM2}
|\hat\cR^\eps\big(f^\eps|\tilde \bw_\eps|^2\big)-\hat\cR^0\big(f^0|\tilde \bw_0|^2\big)|_{-3/2+\kappa;\,T}
\lesssim\vn{f^\eps;f^0}_{\gamma,w;\,T}+o(1)
\end{equ}
holds as $\eps\to 0$, uniformly in modelled distributions $f^\eps\in\cD^{\gamma,w}(V,\Gamma^\eps)$, $f^0\in\cD^{\gamma,w}(V,\Gamma^0)$ bounded by a constant $C$, and in $T\in(0,1]$;
\item For $\eps>0$, $f\in\cD^{\gamma,w}(V)$, one has the identity
\begin{equ}\label{eq:hatR-PAM3}
\hat\cR^\eps\big(f|\tilde \bw_\eps|^2\big)=
\cR^\eps\big(f|\tilde \bw_\eps|^2\big)
%(\cR^\eps f)\big(|\Psi_\eps|^2-C_\eps^{\<PAM-2a-Small>}\big)-C_\eps^{\<PAM-4a-Small>}\scal{f,\<PAM-I2a>}
-\big(a_\rho+\tfrac{|{\log \eps}|}{8\pi}\big)\delta_\d (\cR^\eps f),
\end{equ}
where $a_\rho$ is as in Lemma \ref{lem:reno Psi^2 PAM}.
\end{claim}
\end{lemma}
\begin{proof}
By the definition of the model and Lemma \ref{lem:reno Psi^2 PAM}, for $\eps>0$ one has the identity
\begin{equs}
\cR^\eps\big(f|\tilde \bw_\eps|^2\big)%&=(\cR^\eps f)|\Psi_\eps|^2
%&=(\cR^\eps f)\big(|\Psi_\eps|^2-\E|\Psi_\eps|^2\big)+(\cR^\eps f)\E|\Psi_\eps|^2
&=\cR^\eps(f\<PAM-2new>)+(\cR^\eps f)\delta_\d  R^1_\eps+(\cR^\eps f)\scR R^2_\eps+(\cR^\eps f)R^3_\eps
\\
&\qquad+(\cR^\eps f)\big(a_\rho+\tfrac{|{\log \eps}|}{8\pi}\big)\delta_\partial.\label{eq:decomposing}
\end{equs}
Therefore it suffices to show that the first four terms on the right-hand side extend in an appropriately continuous way to $\eps=0$.

First note that if $f\in\cD^{\gamma,w}(V)$, then $f\<PAM-2new>\in\cD^{\gamma',w'}$ with $\gamma'=\gamma-1-2\kappa=1/2+\kappa$ and $w'=w+(-1-2\kappa)_3$, that is, $\eta'=-3/2+\kappa$, $\sigma'=-5\kappa$, $\mu'=-3/2+\kappa$.
Although one has $\sigma>-1$, unfortunately $\alpha'=-1-2\kappa<1$, so
there is no canonical reconstruction on $\cD^{\gamma',w'}$ as \cite[Sec~4.3]{GH17} does not apply.
This barrier is exactly what is overcome in \cite[App~C]{HP19}:
since in our models we have a decomposition $\Pi_x^\eps=\Pi_x^{+,\eps}+\Pi_x^{-,\eps}$, by \cite[Thm~C.5]{HP19}, there exists maps
$\hat \cR^\eps_0:\cD^{\gamma',w'}\to\CC^{-3/2+\kappa}$ satisfying \eqref{eq:hatR-PAM1} and the desired continuity properties.
Furthermore, for $\eps>0$ one has the identity
$\hat \cR^\eps_0(f\<PAM-2new>)=\cR^\eps(f\<PAM-2new>),$ therefore indeed extending the first term on the right-hand side of \eqref{eq:decomposing}.

For the other terms we may forget that $f$ is a modelled distribution and consider $\cR^\eps f$ as a generic element $F$ of $\CC^{1-2\kappa,w}_P$.
By Proposition \ref{prop:weighted holder2} (a)  $\Tr_{\R_+\times \d D}$ maps $\CC^{1-2\kappa,w}_P$ continuously into
$\CC^{\sigma,\eta\wedge\mu}_{P_0}(\R_+\times\d D)=\CC^{1-3\kappa,-1/2+3\kappa}_{P_0}(\R_+\times\d D)$.
By Proposition \ref{prop:weighted holder2} (b) multiplication on the latter space with $\CC^{1-\kappa,-\kappa}_{\d^2}(Q_\d)$ is continuous in both arguments to $\CC^{1-3\kappa,(-1/2+3\kappa,-\kappa,-1/2+2\kappa)}_{P_0,\d^2}(\R_+\times\d D)$, which in turn is continuously embedded into $\CC^{-1/2+2\kappa}(\R\times\d D)$ by Proposition \ref{prop:weighted holder2} (c).
Finally, by Proposition \ref{prop:weighted holder extend}, $\CC^{-1/2+2\kappa}(\R\times\d D)$ continuously embeds into
$\CC^{-3/2+2\kappa}(\R^4)$. This handles the second term on the right-hand side of \eqref{eq:decomposing}.

Next, we claim that for $F\in\CC^{1-2\kappa,w}_P(\R_+\times D)$, and $R\in\CC_\d^{1-\kappa,-1-\kappa}(Q)$ such that $D_{x_1}R=D_{x_2}R=0$, the distribution $(\scR R) F$ defined as
\begin{equ}
\scal{(\scR R) F,\varphi}=\int_{\R_+\times D}R(z)\Big(F\varphi(z)-F\varphi(\pi_\d z)\Big)\,dz
\end{equ}
is meaningful for all $\eps\in[0,1]$, belongs to $\CC^{-3/2+2\kappa}$,  and is a continuous function of both arguments.
First notice that away from the boundaries this is just a continuous function, which furthermore has an integrable blowup at $P_0$.
In particular, the behavior on the set $\{z:|z|_{P_0}\leq2|z|_{P_1}\}$ is easily controlled as in \cite[Sec.~4.3]{GH17}, and so we focus on test functions centered on the region $\{z:|z|_{P_0}\leq2|z|_{P_1}\}$ and supported away from $P_0$.
%Furthermore, also as in [8], we may restrict our attention to $C_m:=\{z:|z|_{P_1}\leq2|z|_{P_0},|z|_{P_0}\in[2^{-m},2^{-m+2}]\}$, for a fixed (but arbitrary) $m\in\mathbb{N}$.
Take $R$ and $F$ as above with norm $1$, $y\in\{z:|z|_{P_0}\leq2|z|_{P_1}\}$ and a normalised test function $\varphi$.
First consider the case $\lambda\leq |y|_{P_1}/2$.
Then we can simply write
\begin{equ}
|\scal{(\scR R) F,\varphi_y^\lambda}|=|\scal{R F,\varphi_y^\lambda}|
\leq \sup_{z\in\supp\varphi_y^\lambda}|RF(z)|\lesssim |y|^{-1-\kappa}_{P_1}|y|_{P_0}^{-1/2+3\kappa}\leq \lambda^{-3/2+2\kappa},
\end{equ}
which is a bound of required order.
In the case $|y|_{P_1}\leq 2\lambda$ the boundary terms have to be taken into account. One has
\begin{equs}
|\scal{(\scR R) F,\varphi_{y}^\lambda}|&\lesssim
\int_{\supp\varphi_{y}^\lambda}|R(z)|\big|F(z)(\varphi(z)-\varphi(\pi_\d z))+(F(z)-F(\pi_\d z))\varphi(\pi_\d z)\big|\,dz
\\
\big(|z|_{P_0}^{-1/2+3\kappa}\lambda^{-6}|z|_{P_1}+|z|_{P_0}^{-3/2+6\kappa}|z|_{P_1}^{1-3\kappa}\lambda^{-5}\big)\,dz
\\
&\lesssim |y|_{P_0}^{-1/2+3\kappa}\lambda^{-1-\kappa}+|y|_{P_0}^{-3/2+6\kappa}\lambda^{-4\kappa}\lesssim \lambda^{-3/2+2\kappa},
\end{equs}
where we have used that since the exponent of $|z|_{P_1}$ is above $-1$, the integral is finite.
Therefore, $(\scR R)F$ behaves as required.

Finally, for the fourth term, it suffices to notice that $\CC^{1-2\kappa,w}_P$ embeds continuously into $\CC^{-1/2+3\kappa}(\R^4)$,
on which by Proposition \ref{prop:weighted holder mult} multiplication with $\CC^{1-\kappa,-1-\kappa}_{\d^2}(Q)$ is continuous and maps into $\CC^{-1/2,-3/2+2\kappa}_{\d^2}(Q)$, which in turn by Proposition \ref{prop:weighted holder extend} embeds continuously into $\CC^{-3/2+2\kappa}(\R^4)$.
\end{proof}
Using $\hat\cR^\eps$ from Lemma \ref{lem:hatR-PAM} in the operator $\scP_\eps$, the fixed point problem \eqref{eq:PAM-abstract} is completely well-defined and one can use \cite[Thm~5.6]{GH17} to conclude that,
provided $\kappa>0$ is sufficiently small, \eqref{eq:PAM-abstract} has a unique local solution $\bv_\eps$.
Now the main statements follow easily.

\begin{proof}[Proof of Theorem \ref{thm:PAM-N}]
Unlike in the Dirichlet case, the difference between $\hat \cR^\eps$ and the canonical extension of $\cR^\eps$ \emph{does} effect the equation.
Using \eqref{eq:hatR-PAM3}, we have that the
function $v_\eps:=\cR^\eps\bv_\eps$ do not exactly satisfy \eqref{e:etPAM},
but rather
\begin{equ}
(\partial_t -\Delta )v_\eps
=v_\eps(|\Psi_\eps|^2-C_\eps)- 2\nabla v_\eps\cdot\Psi_\eps
-\big(a_\rho+\tfrac{|{\log \eps}|}{8\pi}\big)\delta_\d  v_\eps,
\end{equ}
with $0$ Neumann boundary condition.
That is however equivalent to \eqref{e:etPAM} with boundary condition
\begin{equ}
\partial_n v_\eps=-(a_\rho+\tfrac{|{\log \eps}|}{8\pi})v_\eps.
\end{equ}
Hence, $v_\eps e^{-Y_\eps}$ satisfies \eqref{e:ePAM} with boundary condition
\begin{equ}
\partial_n (v_\eps e^{-Y_\eps})
=-(\partial_n Y_\eps)(v_\eps e^{-Y_\eps})+(\partial_nv_\eps) e^{-Y_\eps}=
-\big(a_\rho+\tfrac{|{\log \eps}|}{8\pi}\big)v_\eps e^{-Y_\eps},
\end{equ}
and so it coincides with $u_\eps$.
The proof is then finished as above.
\end{proof}

\subsection{\texorpdfstring{$\Phi^4_3$}{Phi\^4\_3}}
In the case of the $\Phi^4_3$, one has to keep in mind that the remainder kernel may change as $\eps\to 0$.
As before, we start by defining the modelled distributions,
for $(\eps,a)\in[0,1]\times[0,\infty]$,
\begin{equ}
\bw_{\eps,a}=\scP_{\eps,a}(\bone_D\<0>).
\end{equ}
The reconstruction of $\bone_D\xi_\eps$ is again straightforward. As for $Z_\eps$, in the notation of Lemma \ref{lem:robin}, we choose it to be $Z^{(3a)}$.
As before, it is straightforward that $\bw_{\eps,a}$ belongs to $\cD^{1+3\kappa,(-1/2-\kappa)_3}$, it
depends continuously on $(\eps,a)\in[0,1]\times[0,\infty]$,
and that $\cR^\eps\bw_{\eps,a}=\Psi_{\eps,a}$.

Similarly to Section \ref{sec:restrict}, one can see that $\Psi_{\eps,a}$ can be restricted to the \emph{temporal} hyperplane $P_0=\{t=0\}$
(in fact this is already done in the case without boundaries in [10]).
Let again be $\bi_{\eps,a}$ be the modelled distribution
responsible for the contribution of the initial condition:
it is simply obtained by taking the Taylor polynomial lift of
the function $t,x\mapsto\scal{(\bar\cK+Z^{(3a)})_t(x,\cdot),u_0-\Psi_{\eps,a}\restriction_{P_0}}$.
Equation \eqref{eq:Phi4-lin} then motivates to consider the fixed point problems
\begin{equ}\label{eq:Phi-abstract}
\bv_{\eps,a}=-\scP_{\eps,a}\big(\bone_D^+(\bv_{\eps,a}^3+3\bv_{\eps,a}^2\bw_{\eps,a})\big)
\end{equ}
The operators $\scP_{\eps,a},\scP_{\eps,a}',\scP_{\eps,a}''$ are built from the models $(\Pi^{\eps,a},\Gamma^\eps)$ and the remainder $Z^{(3a)}$, the only difference in them will be the choice of $\hat\cR^\eps$.

We now specify the spaces of modelled distributions on which \eqref{eq:Phi-abstract} can be shown to be a contraction and therefore to have a unique solution.
Denote by $T_0^{\geq \alpha}$ the subspace of $T_0$ generated by basis vectors of homogeneity greater or equal to $\alpha$
(note that no new symbols are included in these spaces).
Fix the exponents $\gamma=1+3\kappa$, $\gamma_0=\kappa$,
$\alpha=0$, $\alpha_0=-1/2-\kappa$, $\alpha_1=-1-2\kappa$, $\alpha_2=-3/2-3\kappa$, the sectors
\begin{equ}
V=\cI(T_0^{\geq \alpha_0})+\bar T,\quad V_0=T_0^{\geq \alpha_0},\quad V_1=T_0^{\geq \alpha_1},\quad V_2=T_0^{\geq\alpha_2},
\end{equ}
as well as
\begin{equs}
%\alpha &=0, &\quad
\eta &=-2/3+\kappa&\quad \sigma &=1/2-4\kappa &\quad \mu &=-2/3+\kappa;\\
%\alpha_0 &=-1-2\kappa, &\quad
\eta_0 &=-2+3\kappa,&\quad \sigma_0 &=-1/2-\kappa,&\quad \mu_0 &=-2+\kappa;\\
%\alpha_1 &=-1/2-3\kappa, &\quad
\eta_1 &=-5/3-\kappa,&\quad \sigma_1 &=-1-2\kappa,&\quad \mu_1 &=-5/3-\kappa;\\
\eta_2 &=-3/2-3\kappa,&\quad \sigma_2 &=-3/2-3\kappa,&\quad\mu_2 &=-3/2-3\kappa.
\end{equs}
It is easy to verify that, by \cite[Lem~5.4]{GH17}, $\bi_{\eps,a}\in\cD^{\gamma,w}(V)$ depends continuously on $(\eps,a)\in[0,1]\times[0,\infty]$, provided $\kappa$ is sufficiently small.
It also follows from \cite[Lem.~4.3]{GH17} that the mappings
\begin{equ}
f\mapsto f^3+3f\bw_\eps,\quad\quad f\mapsto f\bw_\eps^2
\end{equ}
are locally Lipschitz continuous from $\cD^{\gamma,w}(V)$
to $\cD^{\gamma_0,w_0}(V_0)$ and
$\cD^{\gamma_0,w_1}(V_1)$, respectively,
and that $\bw_\eps^3\to\bw_0^3$ in $\scD^{\gamma_2,w_2}(V_2)$.
Furthermore, since $\sigma_0\wedge\alpha_0>-1$,
\cite[Thm~4.9]{GH17} provides a continuous extension $\hat \cR^\eps$ of the reconstruction operator on $\cD^{\gamma_0,w_0}(V_0)$,
which agrees with the canonical extension for $\eps>0$.
This now completes the definition of $\scP_\eps$ in \eqref{eq:Phi-abstract}.

Furthermore, in the case of Dirichlet boundary conditions, that is, $a=\infty$, we are in an even better shape. Indeed, $\scP_{\eps,\infty}'$ and $\scP_{\eps,\infty}''$ are also well-defined, thanks to Lemma \ref{lem:Dirichlet-convolution}.
By \cite[Thm~5.6]{GH17} one can conclude that,
provided $\kappa>0$ is sufficiently small, it has a unique local solution $\bv_{\eps,\infty}$. Now the proof can be easily concluded.
\begin{proof}[Proof of Theorem \ref{thm:Phi4-D}]
Invoking the usual arguments concerning the renormalisation in the bulk (see e.g. \cite[Sec~9]{H0}),
one sees that the function $v_{\eps,\infty}=\cR^\eps\bv_{\eps,\infty}$ solves
\begin{equ}
\end{equ}
with $0$ Dirichlet boundary conditions.
Hence $v_{\eps,\infty} +\Psi_{\eps,\infty}$ solves solves \eqref{e:Phi43-eps},
trivially also with $0$ Dirichlet boundary conditions,
and therefore it coincides with $u_\eps^{\Dir}$.
By \cite[Thm.~5.6]{GH17}, $\bv_{\eps,\infty}\to\bv_{0,\infty}$, the latter
clearly does not depend on $\rho$,
and so the limit $\lim_{\eps\to0}u_\eps^{\Dir}=(\cR^0\bv_{0,\infty})+\Psi_{0,\infty}=:u^{\Dir}$ is also independent of $\rho$, yielding the claim.
\end{proof}

\begin{remark}\label{rem:Dirichlet}
To see in which sense $u^{\Dir}$ vanishes on the boundary, first note that $\cR^0\bv_{0,\infty}$ is a continuous function with $0$ trace.
It is also clear that by definition, $\Psi_{\eps,\infty}$ vanishes on the boundary for $\eps>0$.
By Corollary \ref{cor:restrict1}, $\Psi_{0,\infty}$ has a meaningful trace, which is also $0$.
\end{remark}

Now we move on to the setup of Theorem \ref{thm:Phi4-N}.
First consider the case $b\geq 0$.
Take $b_\eps$ as in the theorem and $c_\eps$ as constructed from $b_\eps$ in Lemma \ref{lem:reno Phi43 square}.
Consider the sequence $\bv_{\eps,c_\eps}$.
Since $c_\eps<\infty$, we are out of the scope of Lemma \ref{lem:Dirichlet-convolution}, so we have to construct $\hat \cR^\eps$ for the operators $\scP_{\eps,c_\eps}'$ and $\scP_{\eps,c_\eps}''$.

The former goes in essentially the same way as in the case of the PAM.
We can summarise the relevant properties of $\hat \cR^\eps$ in the following statement. Its proof follows from the same arguments as that of Lemma \ref{lem:hatR-PAM}
(simply using Lemmas \ref{lem:reno Phi43 square} and \ref{lem:Phi model convergence} in place of Lemmas \ref{lem:reno Psi^2 PAM} and \ref{lem:PAM model convergence}, respectively)
and is therefore omitted.

\begin{lemma}\label{lem:hatR-Phi}
Let $a_\rho$, $b_\eps$, and $c_\eps$ as in Lemma \ref{lem:reno Phi43 square}.
Then for all $\eps\in[0,1]$ there exist maps $\hat \cR^\eps:\scD^{\gamma,w}(V)\bw_{\eps,c_\eps}^2\to\CC^{-1-2\kappa}$ with the following properties.
\begin{claim}
\item For all $\eps\in[0,1]$ and $\psi\in\cC^\infty_0(\R_+\times D)$ one has
\begin{equ}\label{eq:hatR-Phi1}
\scal{\hat\cR^\eps g,\psi}=\scal{\cR^\eps g ,\psi};
\end{equ}
\item The bound
\begin{equ}\label{eq:hatR-Phi2}
\lesssim\vn{f^\eps;f^0}_{\gamma,w;\,T}+o(1)
\end{equ}
holds as $\eps\to 0$, uniformly in modelled distributions $f^\eps\in\cD^{\gamma,w}(V,\Gamma^\eps)$, $f^0\in\cD^{\gamma,w}(V,\Gamma^0)$ bounded by a constant $C$, and in $T\in(0,1]$;
\item For $\eps>0$, $f\in\cD^{\gamma,w}$, one has the identity
\begin{equ}\label{eq:hatR-Phi3}
\hat\cR^\eps(f \bw_{\eps,c_\eps}^2)=
\cR^\eps(f \bw_{\eps,c_\eps}^2)
%(\cR^\eps f)\big(|\Psi_\eps|^2-C_\eps^{\<PAM-2a-Small>}\big)-C_\eps^{\<PAM-4a-Small>}\scal{f,\<PAM-I2a>}
-\big(a_\rho+b_\eps+c_\eps\big)\delta_\d (\cR^\eps f).
\end{equ}
\end{claim}
\end{lemma}

Moving on to $\scP_{\eps,c_\eps}''$, notice that
for each $\eps\in[0,1]$ there is only a single modelled distribution, $\bone_D^+\bw_{\eps,c_\eps}^3$, on which $\hat \cR^\eps$ needs to be constructed.
Also notice that one has for $\eps>0$
\begin{equ}\label{eq:cubic1}
\cR^\eps(\bone_D^+\bw_{\eps,c_\eps}^3)=\bone_D^+\big(\Psi_{\eps,c_\eps}^3-3\ell_\eps({\<Phi-2-Small>})\Psi_{\eps,c_\eps}\big).
\end{equ}
Let us then set, for $\eps\in[0,1]$,
\begin{equ}\label{eq:cubic2}
\hat\cR^\eps(\bone_D^+\bw_{\eps,c_\eps}^3)=\bone_D^+\Pi^{\eps,c_\eps}\<Phi-3new>+
3\bone^+\big((\delta_\d  R^1_\eps)\Psi_{\eps,c_\eps}+(\scR R^2_\eps)\Psi_{\eps,c_\eps}+R^3_\eps\Psi_{\eps,c_\eps}\big).
\end{equ}
First notice that by \eqref{eq:cubic1}, \eqref{eq:cubic2}, the definition of $\Pi^{\eps,c_\eps}\<Phi-3new>$, and Lemma \ref{lem:reno Phi43 square}, we have that for $\eps>0$
\begin{equ}\label{eq:cubic3}
\hat\cR^\eps(\bone_D^+\bw_{\eps,c_\eps}^3)=
\cR^\eps(\bone_D^+\bw_{\eps,c_\eps}^3)
-3\bone^+\big(a_\rho+b_\eps+c_\eps\big)\delta_\d (\Psi_{\eps,c_\eps}).
\end{equ}
In particular, $\hat\cR^\eps(\bone_D^+\bw_{\eps,c_\eps}^3)$ and $\cR^\eps(\bone_D^+\bw_{\eps,c_\eps}^3)$ coincide in the interior of $\R_+\times D$, as required.
It is also clear that $\hat\cR^\eps$ has the required continuity properties:
for the first term on the right-hand side of \eqref{eq:cubic2} this follows from the convergence of the models, while for the second this is precisely the content of Lemma \ref{lem:restrict2} (including that the expressions in \eqref{eq:cubic2} are actually meaningful when $\eps=0$).

%The construction of $\hat \cR^\eps$ on $\scD^{\gamma,w}\bw_\eps^2$ goes just as in Lemma \ref{lem:hatR-PAM}. \mateText{TBC. I think it really is almost the same}

\begin{proof}[Proof of Theorem \ref{thm:Phi4-N}]
It follows from the above that the abstract equation \eqref{eq:Phi-abstract} with $c_\eps$ in place of $a$ admits local solutions $\bv_{\eps,c_\eps}\in\cD^{\gamma,w}$.
Furthermore, by \eqref{eq:hatR-Phi3} and \eqref{eq:cubic3}, $v_{\eps,c_\eps}=\cR^\eps\bv_{\eps,c_\eps}$ satisfy \eqref{eq:Phi4-lin} with boundary conditions $\d_nv_{\eps,c_\eps}=-3c_\eps v_{\eps,c_\eps}$.
The argument leading to \eqref{eq:Phi4-lin} then shows that $\bar u_\eps:=v_{\eps,c_\eps}+\Psi_{\eps,c_\eps}$ satisfy \eqref{e:Phi43-eps} with boundary conditions $\d_n\bar u_\eps=3(a_\rho+b_\eps)\bar u_\eps$ and so it coincides with $u_\eps$.
By Proposition \ref{prop:tiny extension} we have that $\bv_{\eps,c_\eps}\to\bv_{0,b}$, the latter clearly does not depend on $\rho$, so the limit $\lim_{\eps\to0}u_\eps=(\cR^0\bv_{0,b})+\Psi_{0,b}=:u^b$ is also independent of $\rho$.
Finally, one notices that $u^{\infty}$ coincides with $u^{\Dir}$ as defined in the proof of Theorem \ref{thm:Phi4-D}.
\end{proof}

The case $b<0$ requires a very slight modification in order to reduce ourselves to the case $a = 0$ in
\eqref{eq:Psi}.  Defining $\tilde b_\eps=b_\eps-b$, the sequence $(\tilde b_\eps)_{\eps\in(0,1]}$ satisfies \eqref{eq:boundary-reno} with $\tilde b=0$.
Replacing the maps $\hat\cR^\eps$ from Lemma \ref{lem:hatR-Phi} and \eqref{eq:cubic2} by $\tilde\cR^\eps$ defined as
\begin{equs}
\tilde \cR^\eps (f\bw_{\eps,\tilde c_\eps}^2)&:=\hat \cR^\eps (f\bw_{\eps,\tilde c_\eps}^2)+b\delta_\d(\cR^\eps f),
\\
\tilde\cR^\eps(\bone_D^+\bw_{\eps,\tilde c_\eps}^3)&:=\hat\cR^\eps(\bone_D^+\bw_{\eps,\tilde c_\eps}^3)+3\bone^+ b\delta_\d\Psi_{\eps,\tilde c_\eps},
\end{equs}
the proof is concluded just as above.

\appendix
\section{Abstract Schauder estimate -- the Dirichlet case}
We wish to obtain a Schauder estimate in the special case of the Dirichlet heat kernel, but \emph{without} dealing with spatial extensions of $\cR f$, even in the regime $\alpha\wedge\sigma\in (-2,-1)$.
We remark that a similar issue was encountered in \cite[Sec~3]{Cyril19}, but therein the problem could be handled with a more ad hoc workaround.
Let us first formulate what we mean by a Dirichlet kernel.
Simplify notations by setting, in the context of [8], $P_0=\{0\}\times\R^{d-1}$, $P_1=\R\times \partial D$, and
\begin{equ}
\beta=2,\, \frm_0=2, \,\frm_1=1.
\end{equ}
Denote furthermore $Q=[-1,1]\times D$.
\begin{assumption}\label{as:Dirichlet}
Let $K$ (with a corresponding abstract integration map $\cI$) and $Z$ be as in [8].
Suppose furthermore that for all $n\in\N$, $G_n:=K_n+Z_n$ is non-anticipative, symmetric, and vanishes on $(\R^d\setminus P_1)\times P_1$.
\end{assumption}
Note that this implies the following: for any $y\in Q$, $D_2^kG_n(\bar x,y)=D_1^kG_n(y,\bar x)$, and the latter vanishes whenever $\bar x$ belongs to the boundary, so
\begin{equ}\label{eq:Dirichlet impr}
\sup_y |D^k_2G_n(x,y)|\lesssim |x|_{P_1}\sup_{y,z}|D^k_1 \nabla_2 G_n(y,z)|\lesssim |x|_{P_1}2^{n(|k|_\frs+|\frs|-1)}.
\end{equ}
So, in the case $|x|_{P_1}\les 2^{-n}$, we improve the usual bound $|D_2^kG_n|\les 2^{n(|k|_\frs+|\frs|-2)}$ to $|D_2^kG_n|\les |x|_{P_1}2^{n(|k|_\frs+|\frs|-1)}$.

Next, take a partition of unity (essentially the same as in the proof of \cite[Prop.~6.9]{H0}) with the following properties:
\begin{itemize}
\item One has the bounds $|D^k\varphi_{m,i}|\lesssim 2^{m|k|_{\frs}}$, $|\supp\varphi_{m,i}|\lesssim 2^{-m|\frs|}$, as well as the property $\supp\varphi_{m,i}\subset\{y:\,|y|_{P_1}\sim 2^{-m}\}$.
\item
$\sum_{m\in \N}\sum_{i=1}^{\ell(m)}\varphi_{m,i}(y)=1$ for all $y\in Q$, and there exists a constant $C$ such that for all $m,i$, the support of $\varphi_{m,i}$ intersects with $\varphi_{m,j}$ for at most $C$ indices $j$ (in particular, $\ell(m)\les 2^{m(|\frs|-1)}$).
\end{itemize}
Let us set
\begin{equ}
\end{equ}
We shall only use these functions when $2^{-n}\geq|x|_{P_1}$, in which case $G_n^{m,i}$ is nonzero only if $n\leq m$.
Also, for each $n\leq m$ there is a set $A_{n,m}\subset\{1,\ldots,\ell(m)\}$ with cardinality of order $2^{(m-n)(|\frs|-1)}$ such that $G_n^{m,i}=0$ for $i\notin A_{n,m}$.
We can now define our candidate for acting on $G(x,\cdot)$ by distributions $\zeta$ on $Q$:
\begin{equ}\label{eq:G in D}
\zeta(G(x,\cdot)):=\sum_{2^{-n}< |x|_{P_1}}\zeta\big(G_n^{}(x,\cdot)\big)
\sum_{2^{-n}\geq |x|_{P_1}}\sum_{m\in\N}\sum_{i\in A_{n,m}}\zeta\big(G_n^{m,i}(x,\cdot)\big),
\end{equ}
which is meaningful provided the infinite sums (the first one in $n$, the second one in $m$) converge. Note also that if $\zeta$ happens to be a piecewise smooth function up to $P_1$, vanishing for negative times, then the following identity holds on $Q_+=[0,1]\times D$:
\begin{equ}
\zeta(G(x,\cdot))=\int_Q G(x,y)\zeta(y)\,dy.
\end{equ}
The improvement in the Dirichlet case will follow from the following simple property.
It is trivial from \eqref{eq:Dirichlet impr} that $G_n^{m,i}(x,\cdot)$ is $|x|_{P_1}2^{n(|\frs|-1)}2^{-m|\frs|}$ times an approximate Dirac-delta, that is, a function of the form $\psi_y^{2^{-m}}$, with some $y\in Q$ and $\psi\in\cB$ as usual.
However, thanks to the vanishing at the boundary, this can be improved as follows.
\begin{proposition}\label{prop:G vanish}
Suppose $2^{-n}\geq \big(2^{-m}\vee|x|_{P_1}\big)$. Define
\begin{equ}
\hat G_n^{m,i}(x,\cdot)=2^{m(|\frs|+1)}|x|_{P_1}^{-1}2^{-n|\frs|}G_n^{m,i}(x,\cdot).
\end{equ}
Then for any $x\in[-1,1]\times D$, the function $\hat G_n^{m,i}(x,\cdot)$ is of the form $\psi_y^{2^{-m}}$, with $\psi\in\cB$ as usual, and $|y|_{P_1}\geq 2^{-m+1}$.
\end{proposition}
\begin{proof}
It suffices to show that (i) $\supp \hat G_n^{m,i}(x,\cdot)\subset\{y:\,|y|_{P_1}\sim 2^{-m}\}$, (ii) $|\supp\hat G_n^{m,i}(x,\cdot)|\lesssim 2^{-m|\frs|}$, (iii) $|D^k\hat G_n^{m,i}(x,\cdot)|\lesssim 2^{m(|\frs|+|k|_\frs)}$.

Property (i) and (ii) follows easily from the properties of $\varphi_{m,i}$.
For property (iii), first write
\begin{equ}
|D^k\hat G_n^{m,i}(x,y)|\lesssim 2^{m(|\frs|+1)}|x|_{P_1}^{-1}2^{-n|\frs|}\sum_{k_1+k_2=k}\sup_{|z|_{P_1}\sim 2^{-m}}|D^{k_1} G_n(x,z)|\sup_z|D^{k_2}\varphi_{m,i}(z)|
\end{equ}
Since by assumption, is $G_n=0$ on $P_1$, one can write, using \eqref{eq:Dirichlet impr}
\begin{equ}
\sup_{|z|_{P_1}\sim 2^{-m}}|D^{k_1} G_n(x,z)|\lesssim 2^{-m}\sup_{|z|_{P_1}\lesssim 2^{-m}}|\nabla D^{k_1} G_n(x,z)|\les 2^{-m}|x|_{P_1}2^{n(|\frs|+|k_1|_\frs)}
\end{equ}
where $\nabla$ above denotes the spatial gradient. Along with the trivial bound $|D^{k_2}\varphi_{m,i}|\les 2^{m|k_2|_\frs}$, this yields (iii) as claimed.
\end{proof}

To illustrate how this bound can be used, consider $\zeta\in\cC^{\alpha}(Q)$ with $\alpha\in(-2,-1)$. Extending it by $0$ to $\R\times D$ is of course not a problem, since the time boundaries have codimension $2$.
%The intuitive reasoning why $\zeta(G(x,\cdot))$ is meaningful would be \emph{since any two extension $\zeta'$ and $\zeta''$ of $\zeta$ to $\R^d$ differs only by some Dirac-like thing on $\R\times\partial D$, it is cancelled by the vanishing of $G$}. More precisely,
One can write by \eqref{eq:G in D} and Proposition \ref{prop:G vanish}
\begin{equs}[eq:calculation0]
|\zeta(G(x,\cdot)|&\leq \sum_{2^{-n}< |x|_{P_1}}2^{-2n}|\zeta\big(2^{2n}G_n(x,\cdot)\big)|
\\
&\qquad+\sum_{2^{-n}\geq |x|_{P_1}}\sum_{m\in\N}\sum_{i\in A_{n,m}}2^{-m(|\frs|+1)}|x|_{P_1}2^{n|\frs|}|\zeta\big(\hat G_n^{m,i}(x,\cdot)\big)|
\\
&\les \sum_{2^{-n}< |x|_{P_1}}2^{-2n} 2^{-n\alpha}
+\sum_{2^{-n}\geq |x|_{P_1}}|x|_{P_1}2^{n|\frs|}\sum_{m\gtrsim n} 2^{(m-n)(|\frs|-1)}2^{-m(|\frs|+1)}2^{-m\alpha},
\end{equs}
where to get the second inequality, we used that in the first line $\zeta$ is always applied to approximate Dirac-deltas lying completely in $\R\times D$.
Note that the infinite sums converge: the exponent of $2^n$ in the first and of $2^m$ in the second, is $-2-\alpha<0$.
\begin{equ}\label{eq:calculation1}
|\zeta(G(x,\cdot)|\les |x|_{P_1}^{2+\alpha}+\sum_{2^{-n}\geq |x|_{P_1}}2^{-n(2+\alpha)}2^{n}|x|_{P_1}\les |x|_{P_1}^{2+\alpha},
\end{equ}
where in the last inequality we used $1+\alpha<0$. So indeed $\zeta(G(x,\cdot)$ is a well-defined function, and in fact we even see the `improvement by order 2' of the weight at the boundary, since $\zeta(G(x,\cdot)$ vanishes on the boundary with speed $2+\alpha$.

Fixing a $\gamma>0$, we then define
\begin{equ}
\scP f(x)=\cI f(x)+\cJ(x)f(x)+Z\Pi_x f(x)+\bar\cN f(x),
\end{equ}
where we set
\begin{equ}
\bar\cN f(x)=\sum_{|k|_\frs<\gamma+2}\frac{X^k}{k!}\big(\cR f-\Pi_x f(x)\big)\big(D_x^kG(x,\cdot)\big),
\end{equ}
and the last term is understood in the sense of \eqref{eq:G in D}.
We then have the following extension of \cite[Lem~4.12]{GH17}.

\begin{lemma}\label{lem:Dirichlet-convolution}
Suppose that Assumption \ref{as:Dirichlet} holds.
Fix $\gamma > 0$, $w = (\eta,\sigma,\mu)$, and let $V$ be a sector of regularity $\alpha$.
%, and set $a = (\eta\wedge\alpha,\sigma\wedge\alpha,\mu)$.
Assume that these exponents satisfy
\begin{equ}
%\quad \mu>-3
\end{equ}
and take
\begin{equation}\label{eq:exponents2}
\bar \gamma = \gamma + 2,\quad
\bar\eta=(\eta\wedge\alpha)+2, \quad
\bar\sigma=(\sigma\wedge\alpha)+2,\quad
\bar\mu\leq(a_\wedge+2)\wedge 0.
\end{equation}
Suppose furthermore that none of  $\bar \gamma$, $\bar\eta$, $\bar\sigma$, or $\bar\mu$ are integers.
Then, for $f\in\cD_P^{\gamma,w}(V)$, one has
$\scP f\in\cD_P^{\bar \gamma,\bar w}$, where $\bar w=(\bar\eta,\bar\sigma,\bar\mu)$.

Furthermore, if $(\bar\Pi,\bar\Gamma)$ is a second model realising $K$ for $\cI$ and $\bar f\in\cD_P^{\gamma,w}(V,\bar\Gamma)$, then for any $C>0$ the bound
\vn{\scP f;\bar\scP\bar f}_{\bar\gamma,\bar w;\frK}\lesssim\vn{f;\bar f}_{\gamma,w;\bar \frK}+\|\Pi-\bar\Pi\|_{\gamma;\bar\frK}+\|\Gamma-\bar\Gamma\|_{\bar\gamma;\bar\frK}
holds uniformly in models and modelled distributions both satisfying $\vn{f}_{\gamma,w;\bar \frK}+\|\Pi\|_{\gamma;\bar\frK}+\|\Gamma\|_{\bar\gamma;\bar\frK}+\|\zeta\|_{a,\bar \frK}\leq C$, where $\bar\frK$ denotes the $1$-fattening of $\frK$.

Finally, the identity
\begin{equation}\label{eq:reco identity}
\cR\scP f(x)= \cR f(G(x,\cdot))
\end{equation}
holds on $Q_+$, for $f$ vanishing for negative times, where the right-hand side is understood in the sense of \eqref{eq:G in D}.
\end{lemma}

\begin{remark}\label{rem:tilde Rf regularity}
Recall the regularity of $\cR f$, for $f\in\cD_P^{\gamma,w}(V)$ as above: It is a distribution on $Q_+:=[0,1]\times D$ that satisfies the bounds
\begin{equ}\label{eq:cR bound 1}
\cR f(\psi_x^\lambda)\les \lambda^{\eta\wedge\alpha}|x|_{P_1}^{\mu-\eta\wedge\alpha}
\end{equ}
on $Q_+\cap\{|x|_{P_0}\leq 2|x|_{P_1}\}$ for $4\lambda\leq|x|_{P_0}$, and
\begin{equ}\label{eq:cR bound 2}
\cR f(\psi_x^\lambda)\les \lambda^{\sigma\wedge\alpha}|x|_{P_0}^{\mu-\sigma\wedge\alpha}
\end{equ}
on $Q_+\cap\{|x|_{P_1}\leq|x|_{P_0}\}$ for $4\lambda\leq|x|_{P_1}$. It then follows (see e.g. \cite[Prop~2.15]{GH17}) that if $\eta\wedge\alpha>-2$, then $\cR f$ has a unique extension as a distribution on $Q$ such that it agrees with $\cR f$ on $Q_+$, it vanishes on $[-1,0]\times D$, and it satisfies \eqref{eq:cR bound 1} on $Q\cap\{|x|_{P_0}\leq 2|x|_{P_1}\}$ for $2\lambda\leq|x|_{P_1}$.
In the sequel we denote this extension also by $\cR f$.
In particular, $\cR f\in\cC^{\eta\wedge\sigma\wedge\mu\wedge\alpha}(Q)$.
\end{remark}

\begin{proof}
Inspecting the proof of \cite[Lem~4.12]{GH17}, there are two instances where an extension of $\cR f$ is required.
On one hand, when establishing the bound (including making sense of the left-hand side)
\begin{equ}\label{eq:int bound 1}
\sum_{2^{-n}\geq|x|_{P_0}}|\cR f(D^l_1G_n(x,\cdot))|\les |x|_{P_1}^{\bar\mu-|l|_\frs}\left(\frac{|x|_{P_0}}{|x|_{P_1}}\right)^{(\bar\eta-|l|_\frs)\wedge0}
\end{equ}
on the set $|x|_{P_0}\leq|x|_{P_1}$, and the symmetric analogue on the set $|x|_{P_1}\leq|x|_{P_0}$.
On the other hand, the extension of $\tilde \cR f$ is used in the estimate, for $|k+l|_s\geq\bar\gamma$,
\begin{equ}\label{eq:int bound 2}
\sum_{2^{-n}\geq|x|_{P_0}}\|x-y\|^{|k+l|_\frs-\bar\gamma}|\cR f(D_1^{k+l}G_n(\bar y,\cdot))|\les|x|_{P_0}^{\bar\eta-\bar\gamma}|x|_{P_1}^{\bar\mu-\bar\eta}
\end{equ}
whenever $x,y,\bar y$ satisfy $2\|x-\bar y\|\leq 2\|x-y\|\leq|x|_{P_0}\leq|x|_{P_1}$, as well as the symmetric analogue in the case $|x|_{P_1}\leq|x|_{P_0}$.
Our goal is therefore to make sense of these expressions and derive the required bounds without extending $ \cR f$ spatially, making use of Assumption \ref{as:Dirichlet}.
We note that, unlike in [8], the two cases $|x|_{P_0}\lessgtr|x|_{P_1}$ have to be treated slightly differently, since the special property of vanishing on the boundary is only assumed on $P_1$.

We start with \eqref{eq:int bound 1}, when $|x|_{P_0}\leq|x|_{P_1}$, first considering the terms
$|x|_{P_0}\leq 2^{-n}\leq|x|_{P_1}$. By \eqref{eq:cR bound 1}, we get
\begin{equ}
%\sum_{|x|_{P_0}\leq 2^{-n}\leq|x|_{P_1}}
|\cR f(D^l_1G_n(x,\cdot))|
\leq
%\sum_{|x|_{P_0}\leq 2^{-n}\leq|x|_{P_1}}
\leq
\end{equ}
If $\bar\eta-|l|_\frs>0$, then summing over $|x|_{P_0}\leq 2^{-n}\leq|x|_{P_1}$ gives the bound
\begin{equ}
\end{equ}
which is of the required order. For $\bar\eta-|l|_\frs<0$ (equality cannot occur by assumption), the summation yields
\begin{equ}\label{eq:00}
\end{equ}
On the scale $|x|_{P_1}\leq 2^{-n}$, we repeat the calculation in \eqref{eq:calculation0}-\eqref{eq:calculation1}, making use of $\cR f\in\cC^{a_\wedge}(Q)$ to obtain
\begin{equs}\label{eq:01}
\sum_{2^{-n}\geq |x|_{P_1}}&\sum_{m\in\N}\sum_{i\in A_{n,m}}|\cR f(D_1^lG_n^{m,i}(x,\cdot))|
\\
&\les \sum_{2^{-n}\geq |x|_{P_1}}|x|_{P_1}2^{n|\frs|}\sum_{m\gtrsim n}2^{(m-n)(|\frs|-1)}2^{-m(|\frs|+1-|l|_\frs)}2^{-m a_\wedge}\les|x|_{P_1}^{2+a_\wedge-|l|_\frs},
\end{equs}
where in the summation with respect to $n$, we used that $1+a_\wedge-|l|_\frs$ is negative, so the exponent of $2^n$ is positive.
Now suppose $|x|_{P_1}\leq|x|_{P_0}$ and consider the scale $|x|_{P_1}\leq 2^{-n}\leq|x|_{P_0}$.
Using again the argument in \eqref{eq:calculation0}-\eqref{eq:calculation1}, this time making use of the regularity of $\cR f$ given in \eqref{eq:cR bound 2}, we get
\begin{equs}
\sum_{|x|_{P_0}\geq 2^{-n}\geq |x|_{P_1}}&\sum_{m\in\N}\sum_{i\in A_{n,m}}|\cR f(D_1^lG_n^{m,i}(x,\cdot))|
\\
&\les \sum_{|x|_{P_0}\geq 2^{-n}\geq |x|_{P_1}}|x|_{P_1}2^n\sum_{m\gtrsim n}2^{-m(2-|l|_\frs)}2^{-m (\sigma\wedge\alpha)}|x|_{P_0}^{\mu-\sigma\wedge\alpha}
\\
&\les|x|_{P_0}^{\bar\mu-\bar\sigma}|x|_{P_1}\sum_{|x|_{P_0}\geq 2^{-n}\geq |x|_{P_1}}
\end{equs}
If $|l|_\frs+1-\bar\sigma>0$, then this gives a bound of order $|x|_{P_0}^{\bar\mu-\bar\sigma}|x|_{P_1}^{\bar\sigma-|l|_\frs}$, which can be seen to be of the right order just as in \eqref{eq:00}.
If $|l|_\frs+1-\bar\sigma<0$, then the above bound is of order $|x|_{P_0}^{\bar\mu-\bar\sigma-1}|x|_{P_1}\leq|x|_{P_0}^{\bar\mu-\bar\sigma}$, again as required. On the scales $2^{-n}\geq |x|_{P_0}$, we can write similarly to \eqref{eq:01},
\begin{equ}
\sum_{2^{-n}\geq |x|_{P_1}}\sum_{m\in\N}\sum_{i\in A_{n,m}}|\cR f(D_1^lG_n^{m,i}(x,\cdot))|
\les|x|_{P_1}|x|_{P_0}^{1+a_\wedge-|l|_\frs}\leq |x|_{P_0}^{2+a_\wedge-|l|_\frs}.
\end{equ}

To prove \eqref{eq:int bound 2}, first recall that in the situation therein $\|x-y\|\les|x|_{P_i}\sim|y|_{P_i}\sim|\bar y|_{P_i}.$ Therefore by \eqref{eq:int bound 1} one gets
\begin{equs}
\sum_{2^{-n}\geq|x|_{P_0}}&\|x-y\|^{|k+l|_\frs-\bar\gamma}|\cR f(D_1^{k+l}G_n(\bar y,\cdot))|
\\
&\les \|x-y\|^{|k+l|_\frs-\bar\gamma}|x|_{P_1}^{\bar\mu-\bar\eta}|x|_{P_0}^{\bar\eta-|k+l|_\frs}
\les |x|_{P_0}^{\bar\eta-\bar\gamma}|x|_{P_1}^{\bar\mu-\bar\eta},
\end{equs}
where the condition $|k+l|_\frs>\bar\gamma\geq\bar\eta$ was repeatedly used. The analogous estimate in the case $|x|_{P_1}\leq|x|_{P_0}$ in this case follows exactly the same way.
\end{proof}

\endappendix

\begin{thebibliography}{BGHZ19}
\expandafter\ifx\csname url\endcsname\relax
\def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\expandafter\ifx\csname href\endcsname\relax
\def\href#1#2{#2}\fi
\expandafter\ifx\csname burlalt\endcsname\relax
\def\burlalt#1#2{\href{#2}{\texttt{#1}}}\fi

[1]
\textsc{I.~Bailleul}, \textsc{F.~Bernicot}, and \textsc{D.~Frey}.
\newblock Space-times paraproducts for paracontrolled calculus, 3d-{PAM} and
multiplicative burgers equations.
\newblock \emph{Annales scientifiques de l{\textquotesingle}{\'{E}}cole normale
sup{\'{e}}rieure} \textbf{6}, no.~51, (2018), 1399--1456.
\newblock
\burlalt{doi:10.24033/asens.2378}{https://dx.doi.org/10.24033/asens.2378}.

[2]
\textsc{Y.~{Bruned}}, \textsc{F.~{Gabriel}}, \textsc{M.~{Hairer}}, and
\textsc{L.~{Zambotti}}.
\newblock {Geometric stochastic heat equations}.
\newblock \emph{ArXiv e-prints} (2019).
\newblock J. Am. Math. Soc, to appear.
\newblock \burlalt{arXiv:1902.02884}{https://arxiv.org/abs/1902.02884}.

[3]
\textsc{Y.~Bruned}, \textsc{M.~Hairer}, and \textsc{L.~Zambotti}.
\newblock Algebraic renormalisation of regularity structures.
\newblock \emph{Invent. Math.} \textbf{215}, no.~3, (2019), 1039--1156.
\newblock \burlalt{arXiv:1610.08468}{https://arxiv.org/abs/1610.08468}.
\newblock
\burlalt{doi:10.1007/s00222-018-0841-x}{https://dx.doi.org/10.1007/s00222-018-0841-x}.

[4]
\textsc{A.~{Chandra}} and \textsc{M.~{Hairer}}.
\newblock {An analytic BPHZ theorem for regularity structures}.
\newblock \emph{ArXiv e-prints} (2016).
\newblock \burlalt{arXiv:1612.08138}{https://arxiv.org/abs/1612.08138}.

[5]
\textsc{I.~Corwin} and \textsc{H.~Shen}.
\newblock Open {ASEP} in the weakly asymmetric regime.
\newblock \emph{Communications on Pure and Applied Mathematics} \textbf{71},
no.~10, (2018), 2065--2128.
\newblock
\burlalt{doi:10.1002/cpa.21744}{https://dx.doi.org/10.1002/cpa.21744}.

[6]
\textsc{K.~Chouk} and \textsc{W.~van Zuijlen}.
\newblock {Asymptotics of the eigenvalues of the Anderson Hamiltonian with
white noise potential in two dimensions}.
\newblock \emph{The Annals of Probability} \textbf{49}, no.~4(2021).
\newblock
\burlalt{doi:10.1214/20-aop1497}{https://dx.doi.org/10.1214/20-aop1497}.

[7]
\textsc{M.~Gerencs\'{e}r}.
\newblock Nondivergence form quasilinear heat equations driven by space-time
white noise.
\newblock \emph{Ann. Inst. H. Poincar\'{e} Anal. Non Lin\'{e}aire} \textbf{37},
no.~3, (2020), 663--682.
\newblock \burlalt{arXiv:1902.07635}{https://arxiv.org/abs/1902.07635}.
\newblock
\burlalt{doi:10.1016/j.anihpc.2020.01.003}{https://dx.doi.org/10.1016/j.anihpc.2020.01.003}.

[8]
\textsc{M.~Gerencs{\'e}r} and \textsc{M.~Hairer}.
\newblock Singular {SPDE}s in domains with boundaries.
\newblock \emph{Probability Theory and Related Fields} \textbf{173}, no.~3,
(2019), 697--758.
\newblock
\burlalt{doi:10.1007/s00440-018-0841-1}{https://dx.doi.org/10.1007/s00440-018-0841-1}.

[9]
\textsc{P.~Gon{\c c}alves}, \textsc{N.~Perkowski}, and \textsc{M.~Simon}.
\newblock Derivation of the stochastic {Burgers} equation with {Dirichlet}
boundary conditions from the {WASEP}.
\newblock \emph{Annales Henri Lebesgue} \textbf{3}, (2020), 87--167.
\newblock \burlalt{doi:10.5802/ahl.28}{https://dx.doi.org/10.5802/ahl.28}.

[10]
\textsc{M.~Hairer}.
\newblock A theory of regularity structures.
\newblock \emph{Inventiones Mathematicae} \textbf{198}, no.~2, (2014),
\newblock \burlalt{arXiv:1303.5113}{https://arxiv.org/abs/1303.5113}.
\newblock
\burlalt{doi:10.1007/s00222-014-0505-4}{https://dx.doi.org/10.1007/s00222-014-0505-4}.

[11]
\textsc{M.~Hairer} and \textsc{C.~Labb{\'e}}.
\newblock A simple construction of the continuum parabolic {A}nderson model on
\newblock \emph{Electron. Commun. Probab.} \textbf{20}, (2015), 11 pp.
\newblock \burlalt{arXiv:1501.00692}{https://arxiv.org/abs/1501.00692}.
\newblock
\burlalt{doi:10.1214/ECP.v20-4038}{https://dx.doi.org/10.1214/ECP.v20-4038}.

[12]
\textsc{M.~Hairer} and \textsc{K.~Matetski}.
\newblock Discretisations of rough stochastic {PDEs}.
\newblock \emph{The Annals of Probability} \textbf{46}, no.~3(2018).
\newblock
\burlalt{doi:10.1214/17-aop1212}{https://dx.doi.org/10.1214/17-aop1212}.

[13]
\textsc{M.~Hairer} and \textsc{{\'E}.~Pardoux}.
\newblock A {W}ong-{Z}akai theorem for stochastic {P}{D}{E}s.
\newblock \emph{J. Math. Soc. Japan} \textbf{67}, no.~4, (2015), 1551--1604.
\newblock \burlalt{arXiv:1409.3138}{https://arxiv.org/abs/1409.3138}.
\newblock
\burlalt{doi:10.2969/jmsj/06741551}{https://dx.doi.org/10.2969/jmsj/06741551}.

[14]
\textsc{M.~Hairer} and \textsc{{\'E}.~Pardoux}.
\newblock Fluctuations around a homogenised semilinear random {PDE}.
\newblock \emph{Arch. Ration. Mech. Anal.} \textbf{239}, no.~1, (2021),
\newblock \burlalt{arXiv:1911.02865}{https://arxiv.org/abs/1911.02865}.
\newblock
\burlalt{doi:10.1007/s00205-020-01574-8}{https://dx.doi.org/10.1007/s00205-020-01574-8}.

[15]
\textsc{M.~Hairer}, \textsc{M.~Ryser}, and \textsc{H.~Weber}.
\newblock {Triviality of the 2D stochastic Allen-Cahn equation}.
\newblock \emph{Electronic Journal of Probability} \textbf{17}, (2012), 1--14.
\newblock
\burlalt{doi:10.1214/ejp.v17-1731}{https://dx.doi.org/10.1214/ejp.v17-1731}.

[16]
\textsc{M.~Hairer} and \textsc{H.~Shen}.
\newblock A central limit theorem for the {KPZ} equation.
\newblock \emph{Ann. Probab.} \textbf{45}, no.~6B, (2017), 4167--4221.
\newblock \burlalt{arXiv:1507.01237}{https://arxiv.org/abs/1507.01237}.
\newblock
\burlalt{doi:10.1214/16-AOP1162}{https://dx.doi.org/10.1214/16-AOP1162}.

[17]
\textsc{C.~Labb{\'e}}.
\newblock {The continuous Anderson hamiltonian in $d\le 3$}.
\newblock \emph{Journal of Functional Analysis} \textbf{277}, no.~9, (2019),
3187 -- 3235.
\newblock
\burlalt{doi:10.1016/j.jfa.2019.05.027}{https://dx.doi.org/10.1016/j.jfa.2019.05.027}.

[18]
\textsc{S.~Parekh}.
\newblock The {KPZ} limit of {ASEP} with boundary.
\newblock \emph{Communications in Mathematical Physics} \textbf{365}, no.~2,
(2018), 569--649.
\newblock
\burlalt{doi:10.1007/s00220-018-3258-x}{https://dx.doi.org/10.1007/s00220-018-3258-x}.

\end{thebibliography}
\bibliographystyle{Martin}

\end{document}


∎

e1e-mail<EMAIL_ADDRESS>e2e-mail:
<EMAIL_ADDRESS>e3e-mail<EMAIL_ADDRESS>e4e-mail:
<EMAIL_ADDRESS>e5e-mail<EMAIL_ADDRESS>e6e-mail:
<EMAIL_ADDRESS>e7e-mail<EMAIL_ADDRESS>e8e-mail:
<EMAIL_ADDRESS>e9e-mail<EMAIL_ADDRESS>
11institutetext: Instituto Galego de Física de Altas Enerxías (IGFAE),
Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain
22institutetext: European Organization for Nuclear Research (CERN), Geneva,
Switzerland

# Prospects on searches for baryonic Dark Matter produced in $b$-hadron decays
at LHCb

Alexandre Brea Rodríguez e1,addr1  Veronika Chobanova e2,addr1  Xabier Cid
Vidal e3,addr1  Saúl López Soliño e4,addr1  Diego Martínez Santos e5,addr1
Titus Mombächer e6,addr1  Claire Prouvé e7,addr1  Emilio Xosé Rodríguez
Fernández e8,addr1  Carlos Vázquez Sierra e9,addr2

(Published: 2 November 2021)

###### Abstract

A model that can simultaneously explain Dark Matter relic density and the
apparent matter-antimatter imbalance of the universe has been recently
proposed. The model requires $b$-hadron branching fractions to Dark Matter at
the per mille level. The $b$-hadrons decay to a dark sector baryon,
$\psi_{\rm{DS}}$, which has a mass in the region $940~{}\rm{MeV}/c^{2}\leq
m_{\psi_{\rm{DS}}}\leq 4430~{}\rm{MeV}/c^{2}$. In this paper, we discuss the
sensitivity of the LHCb experiment to search for this dark baryon, covering
different types of topology and giving prospects for Runs 3 and 4 of the LHC,
as well as for the proposed Upgrade II. We show that the LHCb experiment can
cover the entire mass range of the hypothetical dark baryon.

††journal: Eur. Phys. J. C

##  1 Introduction

Despite the long standing success of the Standard Model (SM) of particle
physics, there is experimental evidence for New Physics that is not explained
by the SM. The most tantalising signs for New Physics are the apparent
existence of Dark Matter and the matter-antimatter asymmetry in the universe.
Recently, a model has been proposed that can simultaneously explain these two
pressing unknowns Elor:2018twp ; Alonso-Alvarez:2021qfd . In this model,
neutral $B$ mesons oscillate and then decay into a neutral dark-sector baryon
(hereafter called $\psi_{\rm{DS}}$) plus SM hadrons. The matter-antimatter
imbalance is then proportional Elor:2018twp  to

$\displaystyle\begin{split}Y_{B}\propto&\sum_{s,d}A_{SL}^{s,d}\times\mathcal{B}(B^{0}_{s,d}\rightarrow\psi_{\rm{DS}}X)=\\\
&-\sum_{s,d}\bigg{|}\frac{\Delta\Gamma_{s,d}}{\Delta
m_{s,d}\cos(\phi_{s,d}-\arg(\Gamma_{12}^{s,d}))}\bigg{|}\sin(\phi_{s,d}-\arg(\Gamma_{12}^{s,d}))\\\
&\times\mathcal{B}(B^{0}_{s,d}\rightarrow\psi_{\rm{DS}}X)\ ,\end{split}$

where the semileptonic asymmetry is defined as

$A^{s,d}_{SL}=\frac{\Gamma(\overline{B}^{0}_{s,d}\to
f)-\Gamma({B}^{0}_{s,d}\to\bar{f})}{\Gamma(\overline{B}^{0}_{s,d}\to
f)+\Gamma({B}^{0}_{s,d}\to\bar{f})}$

with a final state $f(\bar{f})$ that is specific to
${B}^{0}_{s,d}(\overline{B}^{0}_{s,d})$. It has been directly measured by
several experiments, being the world averages $A^{d}_{SL}=(-21\pm 17)\times
10^{-4}$ and $A^{s}_{SL}=(-6\pm 28)\times 10^{-4}$ HFLAV18 . Alternatively, it
can be indirectly inferred through the $B^{0}_{s,d}$ mixing parameters $\Delta
m_{s,d}$, $\Delta\Gamma_{s,d}$, $\phi_{s,d}$ and $\Gamma^{s,d}_{12}$. In the
SM as well as in New Physics models with negligible contributions to $\Delta
F=1$ penguins, these asymmetries are predicted to be $A^{d}_{SL}=(-4.73\pm
0.42)\times 10^{-4}$ and $A^{s}_{SL}=(0.205\pm 0.018)\times 10^{-4}$
Lenz:2019lvd , where similar values are found in Ref. Asatrian:2020zxa . Note
that, since one needs a positive value for the semileptonic asymmetry in order
to obtain a matter-dominated universe, the possibility of having the
baryogenesis mechanism lead by $B_{s}^{0}$ decays to dark-sector baryons is
significantly favoured over those of $B^{0}$ mesons.

In order to fulfill the observed matter-antimatter imbalance, the inclusive
branching fractions of $b$ hadrons decaying into the dark-sector baryon must
be relatively large, $\gtrsim 10^{-4}$, with experimental upper limits ranging
from $\leq 10^{-4}$ to $\leq 10^{-2}$ depending on the considered
$\psi_{\rm{DS}}$ mass Alonso-Alvarez:2021qfd ; ALEPH:2000vvi . Relatively high
values could in principle also modify $\Gamma_{s}/\Gamma_{d}$ if they have
different values for $B_{s}^{0}$ mesons compared to $B^{0}$ mesons111G.
Alonso, private communication..

As for the exclusive branching fractions (as are the examples shown in this
paper), these are expected to be at the $\gtrsim 10^{-6}$ level. Since the
$\psi_{\rm{DS}}$ baryon must be produced on shell in $B$-meson decays together
with at least one proton (to preserve baryon number), its mass cannot be
larger than $m_{B_{s}^{0}}-m_{p}$. In addition, to prevent the proton decay,
the $\psi_{\rm{DS}}$ baryon must be heavier than the proton. Hence its mass
range is limited to $940\,\text{MeV}/c^{2}\leq m_{\psi_{\rm{DS}}}\,\leq
4430~{}\text{MeV}/c^{2}$. As an additional remark, we point the reader to
Refs. Elor:2020tkc ; Elahi:2021jia , which approach the same problem from
slightly different angles.

Due to its unique design and the abundant production of $b$ hadrons in high-
energetic $pp$ collisions, the LHCb experiment is ideally suited for searching
for the proposed type of baryonic dark matter in $b$-hadron decays
Chobanova:2020vmx ; Borsato:2021aum . Equipped with a highly granular vertex
detector very close to the $pp$-interaction point LHCbCollaboration:2013bkh ,
a precise reconstruction of the $b$-hadron decay vertex can be achieved, which
is crucial to obtain a good sensitivity in the search for the dark baryon. Due
to baryon-number conservation, searches using $B$ mesons need at least one
proton and at least another charged particle to reconstruct the decay vertex.
This implies $B$ mesons alone cannot cover the entire phase space of the
model. However, LHCb can profit from its large production of $b$ baryons,
which do not require the presence of a proton in the final state to preserve
baryon number, and which are also heavier than $B_{s}^{0}$ mesons.

In this paper, we use as benchmark modes the decays
$B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)(\rightarrow pK^{-})$,
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}(\rightarrow\pi^{+}\pi^{-}\Lambda_{c}^{+}(\rightarrow
pK^{-}\pi^{+}))$, $\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}K^{+}\pi^{-}$, and
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}\pi^{+}\pi^{-}$222Charge-conjugate
modes are included unless explicitly stated. for the search of the dark baryon
$\psi_{\rm{DS}}$ with the LHCb experiment. These channels are sensitive to
different $\mathcal{O}_{q_{1}q_{2}}$ operators out of those introduced in
Refs. Elor:2018twp ; Alonso-Alvarez:2021qfd . More precisely,
$B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$ decays would be sensitive to
$\mathcal{O}_{us}$, $B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ to
$\mathcal{O}_{cd}$, $\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}K^{+}\pi^{-}$ to
$\mathcal{O}_{us}$, and
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}\pi^{+}\pi^{-}$ to
$\mathcal{O}_{ud}$. Since none of these operators are favored a priori by the
mechanism, it becomes essential to probe as many of them as possible.

We explore the possibility of further background suppression by tagging some
of these decays via the decay chains
$\Sigma_{b}^{(*)\pm}\rightarrow\Lambda_{b}^{0}\pi^{\pm}$. The projections are
based on the detector geometry proposed for the data taking periods Run 3 and
Run 4 of the LHCb experiment Bediaga:2012uyd  and its Upgrade II Aaij:2636441
. The expected sensitivities are computed at corresponding benchmark
luminosities of $15$ (Run 3), $50$ (Run 4), and $300\,\rm{fb}^{-1}$ (Upgrade
II) Aaij:2636441  of $pp$ collisions at a center-of-mass energy of
$14\,\rm{TeV}$.
The document is organised as follows: in Sec. 2 we describe the tools used to
simulate events and specific decays within the LHCb detector; Section 3
discusses potential event selections; and in Sec. 4 we describe the estimated
sensitivities for the different modes; finally, potential systematic effects
are discussed in Sec. 5; we conclude in Section 6.

##  2 Event generation and simulation

Proton-proton collisions are simulated with Pythia 8.226 pythia  at a center-
of-mass energy of 14 TeV. Note that this does not include $B^{+}_{c}$ mesons
and therefore backgrounds from $B_{c}^{+}$ mesons are not considered in this
study.

The four-momenta and origin point of the particles of interest that have been
produced by Pythia are stored for further processing. Subsequently, signal $b$
hadrons are decayed according to phase space. No bremsstrahlung is included in
the generation of the decays. Background events are obtained from $b\bar{b}$
events in Pythia filtered with selection requirements before proceeding with
detector simulation, in order to reduce the amount of events to process.

To obtain the background and signal yields, we use the $b\bar{b}$ cross
sections as measured by the LHCb experiment at $\sqrt{s}=13$ TeV and
$\sqrt{s}=7$ TeV Aaij:2016avz  extrapolated linearly with the center-of-mass
energy to $\sqrt{s}=14$ TeV. Hadronisation fractions for $b$ hadrons are
obtained from Refs. Aaij:2021nyr  and Aaij:2016avz , assuming $f_{u}=f_{d}$.
The production of heavier $b$ baryons and $b$ hadrons with multiple heavy
quarks is neglected. We also neglect the effect of pile-up, given the dominant
background is composed of single $b\bar{b}$ events and also that we expect the
effect of wrong primary vertex association333i.e. to associate the signal
candidate to a wrong primary vertex, which potentially affects the
determination of variables such the impact parameter. to barely affect the
sensitivity, given the future vertexing detector of LHCb will be 4D
Aaij:2244311  (including timing), which is known to be very useful to mitigate
this effect.

The detector is simulated using the code in Ref. Chobanova:2020vmx . The
simulated detector elements are the RF foil, the VeloPix (VP) stations, the
Upstream Tracker, the Magnet, and the SciFi Tracker LHCbCollaboration:2013bkh
; LHCbCollaboration:2014tuj , and are implemented in python2. No calorimetry
or particle identification (PID) is simulated at this stage. The particles
generated by the procedure described above are passed through the detector
elements and yield hits where appropriate. The set of hits is then processed
by a track fit algorithm which calculates the track slopes at origin, the
momentum, and a point in the early stage of the particle trajectory. Note that
through this procedure all detector acceptance effects are accounted for. The
simulation neglects occupancy effects and hit inefficiency. To account for
this, when obtaining absolute efficiencies, the tracking efficiency is
artificially scaled by $98\%$. The detector simulation reproduces at first
order the impact parameter and momentum resolutions of the existing full
simulation of the LHCb Upgrade, as well as their kinematic dependencies
LHCbCollaboration:2014tuj ; LHCbCollaboration:2013bkh .

The final computation of the expected exclusions for the decay channels under
evaluation is performed by means of the “missing” transverse momentum
($p_{\rm{T}}^{\rm{miss}}$). This is defined as the sum of all the momenta of
the reconstructed daughters in the direction transverse to the $b$-hadron
direction of flight, as illustrated in Fig. 1.

Figure 1: Determination of $p_{\rm{T}}^{\rm{miss}}$ in
$X_{b}\rightarrow\psi_{\rm{DS}}h_{1}^{SM}h_{2}^{SM}$ decays. For the
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ channel, where more
objects are present in the final state, the computation would be analogous. In
the figure, PV (SV) represent the primary (secondary) vertices, i.e., the
$b$-hadron production (decay) points.

This quantity can only be determined if both the $b$-hadron origin and decay
position are known. Both of these properties can be computed experimentally,
since they correspond to the $pp$-interaction point, or the origin of all the
final-state charged tracks of the $b$ hadron, respectively. The requirement to
know the $b$-hadron decay position limits the amount of channels one can
study, since this involves knowing the trajectory of the final state SM
daughters, which is not possible if these are stable or very long-lived
neutral particles. The $p_{\rm{T}}^{\rm{miss}}$ variable is used in our study
to discriminate signal from background. For the
$B\rightarrow\psi_{\rm{DS}}\Lambda$ analyses, the actual expected limit is
determined based on this quantity. Furthermore, for the
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}h^{+}_{1}h^{-}_{2}$ analyses, also
the $h^{+}_{1}h^{-}_{2}$ invariant mass, $m(h^{+}_{1}h^{-}_{2})$, is very
useful since it shows sharp kinematic end points at the $\psi_{\rm{DS}}$ mass
Chobanova:2020vmx , so the limit determination is based on a region selected
in the two-dimensional $(m(h^{+}_{1}h^{-}_{2})-p_{\rm{T}}^{\rm{miss}})$ plane.

##  3 Event selection

All tracks involved in the analysis are required to satisfy
$p_{\rm{T}}>800~{}\textrm{MeV}/c$ and $0.1<\textrm{IP}<3$ mm, where
$p_{\rm{T}}$ is the transverse momentum in the laboratory frame and IP the
impact parameter with respect to the $pp$ interaction vertex444The lower IP
requirement discriminates against particles produced at the $pp$ interaction
vertex, while the IP upper requirement helps to remove potential ghost
tracks.. The only exception to this are the pions appearing in the
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ decay for which the
$p_{\rm{T}}$ requirement is relaxed to 250 MeV$/c$. We assume 100 % trigger
efficiency. The rest of the requirements are channel dependent and are
explained in the next subsections.

###  3.1 Isolation

The main background for the targeted decay channels in this analysis is
composed of hadrons coming from a $b\bar{b}$ pair. A key feature of this
situation is that the background tracks are typically accompanied by other
objects, while for the signal the final state is “isolated” and no other
nearby particles are expected. Based on this feature, we build a quantity that
provides discrimination between signal and background. For this, we take all
charged particles arising from $b\bar{b}$ pairs (having excluded those forming
the signal candidates), and we select only those reconstructed using our fast
simulation procedure, with $p_{\rm{T}}>250$ MeV$/c$. With this, and following
Refs. CidVidal:2019qub ; CidVidal:2019urm , we determine the smallest distance
of closest approach (DOCA) to our candidates and use that to measure the
isolation. The distribution of this quantity, $\mathrm{DOCA_{ISO}}$ will peak
at smaller quantities for background than for signal. The requirement we
apply, $\mathrm{DOCA_{ISO}}>0.05$, retains $\sim 80\%$ of signal candidates
and suppresses $\sim 75\%$ of the background.

###  3.2 Selection for $B\rightarrow\psi_{\rm{DS}}\Lambda$ decays

The selection of each $B\rightarrow\psi_{\rm{DS}}\Lambda$ decay modes is based
on the full reconstruction of the $\Lambda$ resonances, applying invariant
mass and DOCA requirements to the relevant daughter particles. Furthermore, IP
and $p_{\rm T}$ requirements are applied to the $\Lambda$ candidates. The full
list of requirements for these two channels can be found in Table 1, and the
main considerations for each of the selections in the next paragraphs.

For the $B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$ decay mode, the
$\Lambda(1520)$ resonance is reconstructed through its decay to a $pK^{-}$
pair. Each member of the pair is required to be close to each other (applying
a DOCA requirement), and an invariant mass requirement to the pair is also
applied, so the mother candidate has an invariant mass consistent with that of
the $\Lambda(1520)$. Furthermore, the $\Lambda(1520)$ candidate must not be
consistent with originating at the $pp$ interaction vertex (IP requirement)
and satisfy a minimum $p_{\rm T}$ requirement.

For the $B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ channel, two
subsequent decays are required to reconstruct the signal candidates,
$\Lambda_{c}(2595)^{+}\rightarrow\Lambda_{c}^{+}\pi^{+}\pi^{-}$ and
$\Lambda_{c}^{+}\rightarrow pK^{-}\pi^{+}$. After applying the track-level
requirements described above, the $\Lambda^{+}_{c}$ candidates are assumed to
be very clean and no further requirements are applied to them. No fake
combinations of $\Lambda^{+}_{c}$ are considered either, since this is a very
clean resonance at LHCb, being narrow in invariant mass, displaced555In this
case there is a “double” displacement given these $\Lambda^{+}_{c}$ are
produced in $B$-meson decays. and with a final state including a proton, a
kaon and a pion, which are separable from each other through PID requirements.
Real reconstructed $\Lambda^{+}_{c}$ baryons are then combined with pions to
form $\Lambda_{c}(2595)^{+}$ candidates. Then, as in the $\Lambda(1520)$ case,
invariant mass, DOCA, IP and $p_{\rm T}$ requirements are applied to these.

One relevant exclusive background we have considered for both channels is that
originating from $X_{b}\rightarrow Y_{\rm SM}\Lambda_{c}(2595)^{+}$ and
$X_{b}\rightarrow Y_{\rm SM}\Lambda(1520)$ decays, i.e., cases in which the
baryons accompanying $\psi_{\rm{DS}}$ are actually produced in $b$-hadron
($X_{b}$) SM decays together with other SM particles ($Y_{\rm SM}$). Although
closer to signal than the combinatorial background, these type of candidates
are still separable by means of the isolation and $p_{\rm{T}}^{\rm{miss}}$.
After all selection requirements, they are seen to contribute $<0.5\%$
compared to the combinatorial background for the
$B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$ channel, so are not considered
further in this case. On the contrary, for the
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ channel, they are
relevant and become the dominant source of background, so they are accounted
for in the sensitivity determination. The main contribution to this type of
background in this channel is composed of semileptonic $\Lambda_{b}$ decays,
i.e $\Lambda_{b}\rightarrow\Lambda_{c}(2595)^{\pm}l^{\mp}\nu_{l}$, with
$l^{\mp}$ and $\nu_{l}$ being a charged lepton and corresponding neutrino,
respectively.

Table 1: Requirements for each $B\rightarrow\psi_{\rm{DS}}\Lambda$ decay mode.
These values have been chosen in order to maximize the discrimination between
signal and background. Here, “max(DOCA)” refers to the maximum value of all
the DOCA computed for the two-track combinations among the $\Lambda(1520)$ and
$\Lambda_{c}(2595)^{+}$ daughters, and $|m_{\Lambda}-m_{\rm PDG}|$ describes
the mass window in the relevant reconstructed invariant $\Lambda$ mass around
its known mass PDG . The IP and $p_{\rm{T}}$ requirements are applied to the
$\Lambda(1520)$ and $\Lambda_{c}(2595)^{+}$ candidates. All the track
requirements, introduced in Sec. 2, have also been imposed for every decay
mode.

Variable | $B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ | $B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$
---|---|---
$\mathrm{DOCA_{ISO}}$ | 0.05 mm | 0.05 mm
max(DOCA) | 0.15 mm | 0.075 mm
$|m_{\Lambda}-m_{\rm PDG}|$ | 20 MeV/$c^{2}$ | 7.5 MeV/$c^{2}$
IP | 0.07 mm | 0.1 mm
$p_{\rm{T}}$ | 4000 MeV/$c$ | 1000 MeV/$c$

###  3.3 Selection and Multivariate Classifier for
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}h_{1}^{+}h_{2}^{-}$ decays

The $\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}h_{1}^{+}h_{2}^{-}$ decays suffer
similarly from a large amount of combinatorial background from
$b\bar{b}\rightarrow h_{1}^{+}h_{2}^{-}X$ processes. In order to improve the
discrimination power of the analysis, multivariate classifiers are built using
the $p_{\rm{T}}^{\rm{miss}}$ and the $h_{1}^{+}h_{2}^{-}$ reconstructed mass
as input variables. Before the training, the following requirements are
applied: the two meson tracks must satisfy the requirements introduced at the
end of Sec. 2, the DOCA between the two mesons is required to be less than 0.1
mm, the distance of flight in the z direction ($\Delta z$) larger than 5 mm,
the transverse momentum of the track larger than 1 GeV/$c$, and the ratio
between the impact parameter of the meson pair and the distance of flight in
the z direction (IP$(h_{1}^{+}h_{2}^{-})/\Delta z$) is required to be smaller
than 0.1. The classifier relies on the mathematical method used in Ref.
Aaij:2011rja  (see also Ref. doi:10.1063/1.168691 ). A requirement is
performed in the multivariate classifier response such that it minimizes the
expected limit at a given integrated luminosity and $\psi_{\rm{DS}}$ mass.

###  3.4 Efficiencies

The efficiencies of the selection requirements proposed in this paper are
shown in Table 2 for different benchmark masses of the dark baryon,
$m_{\psi_{\rm{DS}}}$.

Table 2: Reconstruction and selection efficiencies for the signal decays
described in the text, as obtained from fast simulation, where we added a
posteriori a $98\%$ efficiency per track to account for multiplicity effects.
Geometry of the detector is considered within $\varepsilon^{REC}$. The
efficiencies are shown in %. The $\varepsilon^{REC\&PT}$ efficiency implies
that all tracks have been required to have a transverse momentum greater than
$800~{}\text{MeV}/c$, except for
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ decays, where this
requirement is relaxed to $250~{}\text{MeV}/c$ for pions. The
$\varepsilon^{SEL/REC}$ efficiency includes multivariate and isolation
requirements when applicable. $\varepsilon^{PID/SEL}$ refers to the efficiency
assumed for the PID requirements, as explained in the text. Note this does not
apply to modes with no protons or kaons in the final state. Finally,
$\varepsilon^{TOTAL}$ is the product of $\varepsilon^{REC}$,
$\varepsilon^{SEL/REC}$ and $\varepsilon^{PID/SEL}$. The number next to
$\psi_{\rm{DS}}$ refers to $m_{\psi_{\rm{DS}}}$, given in MeV$/c^{2}$.

Decay mode | $\varepsilon^{REC}$ | $\varepsilon^{REC\&PT}$ | $\varepsilon^{SEL/REC}$ | $\varepsilon^{PID/SEL}$ | $\varepsilon^{TOTAL}$
---|---|---|---|---|---
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(940)K^{+}\pi^{-}$ | 7.6 | 5.1 | 13.2 | 74.2 | 0.74
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(940)\pi^{+}\pi^{-}$ | 7.3 | 4.8 | 14.3 | - | 1.04
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(1500)K^{+}\pi^{-}$ | 7.7 | 4.8 | 10.8 | 76.0 | 0.63
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(1500)\pi^{+}\pi^{-}$ | 7.3 | 4.5 | 12.5 | - | 0.91
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(2000)K^{+}\pi^{-}$ | 7.7 | 4.5 | 6.67 | 79.2 | 0.41
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(2000)\pi^{+}\pi^{-}$ | 7.4 | 4.1 | 10.9 | - | 0.81
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(2400)K^{+}\pi^{-}$ | 7.8 | 4.1 | 9.11 | 80.8 | 0.57
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(2400)\pi^{+}\pi^{-}$ | 7.4 | 3.7 | 8.89 | - | 0.66
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(4340)K^{+}\pi^{-}$ | 8.1 | 1.2 | 3.12 | 88.0 | 0.22
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(4340)\pi^{+}\pi^{-}$ | 7.4 | 0.97 | 2.88 | - | 0.21
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(4470)K^{+}\pi^{-}$ | 8.2 | 0.91 | 2.18 | 87.7 | 0.16
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}(4470)\pi^{+}\pi^{-}$ | 7.4 | 0.79 | 2.23 | - | 0.17
$B^{+}\rightarrow\psi_{\rm{DS}}(940)\Lambda_{c}(2595)^{+}$ | 5.7 | 1.3 | 20.9 | 56.7 | 0.68
$B^{+}\rightarrow\psi_{\rm{DS}}(1500)\Lambda_{c}(2595)^{+}$ | 5.6 | 1.1 | 19.3 | 56.6 | 0.61
$B^{+}\rightarrow\psi_{\rm{DS}}(2000)\Lambda_{c}(2595)^{+}$ | 5.3 | 0.88 | 18.5 | 57.2 | 0.56
$B^{+}\rightarrow\psi_{\rm{DS}}(2400)\Lambda_{c}(2595)^{+}$ | 4.9 | 0.66 | 16.4 | 57.1 | 0.46
$B^{0}\rightarrow\psi_{\rm{DS}}(940)\Lambda(1520)$ | 13.3 | 12.7 | 38.1 | 56.9 | 2.88
$B^{0}\rightarrow\psi_{\rm{DS}}(1500)\Lambda(1520)$ | 13.7 | 12.9 | 36.6 | 56.4 | 2.83
$B^{0}\rightarrow\psi_{\rm{DS}}(2000)\Lambda(1520)$ | 13.5 | 12.6 | 35.5 | 56.8 | 2.72
$B^{0}\rightarrow\psi_{\rm{DS}}(2400)\Lambda(1520)$ | 13.3 | 6.8 | 34.5 | 56.6 | 2.6
$B^{0}\rightarrow\psi_{\rm{DS}}(3500)\Lambda(1520)$ | 12.1 | 3.5 | 23.6 | 56.5 | 1.61

We also simulate the effect of applying PID requirements that would make the
background contribution from misidentified particles negligible. For kaons, we
take the expected kaon identification efficiency for a Delta Log Likelihood
(DLL) $\rm{DLL}_{K\pi}>5$ requirement as a function of momentum
Collaboration:1624074 , and apply it to our samples, which provides a penalty
efficiency factor. For protons we follow a similar procedure, applying
$\rm{DLL}>5$ requirements to discriminate protons against kaons and pions. In
this case, since no projections exist for the LHCb upgrades, we take the
measured numbers from the data-taking periods Run 1 and Run 2 of the LHCb
experiment Adinolfi:2012qfa .

###  3.5 Tagging

A useful way to deal with background in $\Lambda_{b}^{0}$ decays while getting
additional kinematic information, is to select those $\Lambda_{b}^{0}$
particles coming from $\Sigma_{b}^{(*)\pm}$ as suggested in Ref. Stone:2014mza
. About one third of $\Lambda_{b}^{0}$ baryons are coming from
$\Sigma_{b}^{(*)\pm}$ according to Pythia pythia . More importantly, from Ref.
Aaij:2018tnn  one can infer that LHCb yields approximately one tagged
$\Lambda_{b}^{0}$ per 10 untagged $\Lambda_{b}^{0}$, including the effect of
reconstructing and selecting the associated slow pion from the decay
$\Sigma_{b}^{(*)\pm}\rightarrow\Lambda_{b}^{0}\pi^{\pm}$. Most background
events will have plenty of slow prompt pions coming from the same collision
point as the $\Lambda_{b}^{0}$ candidate. However, background events are not
expected to peak at the $\Sigma_{b}^{(*)}$ mass, while signal events do (see
Fig. 2). This could allow extra selection requirements for further background
reduction, as well as a good signal confirmation due to the distinctive two-
peak pattern. In these sensitivity studies the usage of tagging is not
included, though.

As for the $B\rightarrow\psi_{\rm{DS}}\Lambda$ modes, tagging would also be
possible in principle for the
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ channel through $B^{+}$
mesons produced in $B_{s2}^{\ast 0}\rightarrow B^{+}K^{-}$ decays. Around 20%
of all $B^{+}$ mesons are produced in a $B_{s2}^{\ast 0}$ decay, according to
Pythia. Note that this type of tagging has already been used at LHCb to search
for the lepton flavor violating $B^{+}\rightarrow K^{+}\mu^{-}\tau^{+}$ decay
Aaij:2020mqb , where the $\tau$ lepton was effectively treated as a missing
particle. This search achieved upper limits on the branching fraction of this
decay at the level of $10^{-5}$, using an LHCb data set accounting for an
integrated luminosity of 9 fb-1.

Figure 2: Normalised distributions of the reconstructed invariant mass of
$\Sigma_{b}^{(*)\pm}\rightarrow\Lambda_{b}^{0}(\rightarrow\psi_{\rm{DS}}h^{+}_{1}h^{-}_{2})\pi^{\pm}$
candidates in the simulated samples. The mass of the $\Lambda_{b}$ baryon is a
fixed parameter in this case. Hatched green: signal, blue: $b\bar{b}$
background. Top: $\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}K^{+}\pi^{-}$
candidates. Bottom: $\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}\pi^{+}\pi^{-}$
candidates. The tagging of $\Lambda_{b}^{0}$ allows a clear distinction
between signal and background samples.

##  4 Estimated statistical sensitivities

As discussed in Sec. 2, we obtain the sensitivities from the distributions of
$p_{\rm{T}}^{\rm{miss}}$, or $p_{\rm{T}}^{\rm{miss}}$ versus
$m(h^{+}_{1}h^{-}_{2})$ for signal and background candidates surviving the
selections introduced in Sec. 3. In order to maximize our sensitivity, the
number of observed events in background of the $B$-meson decay modes are
obtained in bins of $p_{\rm{T}}^{\rm{miss}}$, while for the $\Lambda_{b}^{0}$
decay modes a single two-dimensional region as a function of
$p_{\rm{T}}^{\rm{miss}}$ and $m(h^{+}_{1}h^{-}_{2})$ is considered. Figures 3
and 4 show the background and signal distributions (for different values of
the dark baryon mass) for the $B\rightarrow\psi_{\rm{DS}}\Lambda$ and
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}h_{1}h_{2}$ decay modes,
respectively.

Figure 3: Distributions of $p_{\rm{T}}^{\rm{miss}}$ for
$B\rightarrow\psi_{\rm{DS}}\Lambda$ and associated backgrounds. Top:
$B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$. Bottom:
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$. Blue: background from
inclusive $b\bar{b}$ events. Orange: signal for a high mass (close to phase
space limit) dark baryon. Green: signal for a dark baryon at 940 MeV$/c^{2}$.
The signal and background yields correspond to those expected at LHCb at 15
fb-1, with the signal branching fractions assumed to be those that could be
excluded at 95% CL (see Sec. 4). No systematic uncertainties are assumed in
these plots.

Figure 4: Reconstructed invariant mass versus $p_{\rm{T}}^{\rm{miss}}$ for
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}h^{+}_{1}h^{-}_{2}$ candidates. Top:
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}K^{+}\pi^{-}$. Bottom:
$\Lambda_{b}^{0}\rightarrow\psi_{\rm{DS}}\pi^{+}\pi^{-}$. Blue dots:
background from inclusive $b\bar{b}$ events, green circles: signal for a dark
baryon at 940 MeV$/c^{2}$, orange squares: signal for a dark baryon at 4470
MeV$/c^{2}$. Background from light mass resonances, such as
$K^{*}(892)\rightarrow K^{+}\pi^{-}$, $D^{0}\rightarrow K^{+}\pi^{-}$, and
$K_{S}^{0}\rightarrow\pi^{+}\pi^{-}$ is clearly visible in the plots. No
systematic uncertainties are assumed in these plots.

Exclusion limits on the decay branching fractions of the four $B$ and
$\Lambda_{b}^{0}$ decay modes are then computed using these distributions,
with the Modified Frequentist confidence level (CLs) method, described in
Refs. CLs ; CLsv2 . These values are obtained per decay mode and benchmark
$\psi_{\rm{DS}}$ mass hypotheses as listed in Table 2.

Finally, we consider the limits at 95% C.L. per $\psi_{\rm{DS}}$ mass
hypothesis and we interpolate these values to cover intermediate mass points,
as presented in Fig. 5. For $B$-meson decay modes an exponential extrapolation
is used for masses above 3.5(2.5) GeV/c2, to account for the inability to
reconstruct the $\Lambda$ mesons due to the limited phase space.

For all the decay modes, branching fractions down to $(1$ –
$5){\times}10^{-6}$ for a dark baryon mass between 1 and 2.5 GeV/c2, could be
excluded. The best limits are achieved in the low $m_{\psi_{\rm{DS}}}$ region
and are given by the $B^{0}$ mode we used, because of the complex final state.
High masses are only accessible with $\Lambda_{b}^{0}$ decays. It should be
noted that at high masses the differences between the inclusive and exclusive
branching fraction are smaller, and hence higher branching fractions are
expected.

Figure 5: Projected statistical sensitivity in terms of branching fractions
excluded at 95% CL at, from top to bottom, $15~{}\textrm{fb}^{-1}$,
$50~{}\textrm{fb}^{-1}$ and $300~{}\textrm{fb}^{-1}$. All three curves were
obtained by interpolation, having used an exponential extrapolation (dashed
line) in the abrupt change of slope for the $B$ decays in order to account for
the detector’s inability to reconstruct neither the $\Lambda(1520)$ nor the
$\Lambda_{c}(2595)^{+}$ baryons at the limit of the phase space, i.e.,
whenever $m_{\psi_{\rm{DS}}}$ reaches the mass difference of the SM mother
meson and daughter baryon.

##  5 Systematic uncertainties

It should be noted that in these analyses, the background yields are typically
larger than in other searches, such as rare decays. Hence, systematic
uncertainties related to background yield or modelling can potentially limit
the search. For instance, with a background yield systematic uncertainty of
$1\%$ all the searches shown in this study would be already systematically
limited with $1\,\rm{fb}^{-1}$, and with a systematic uncertainty of 1 per
mile, they would be systematically limited after $\sim 15\,\rm{fb}^{-1}$.
Hence, further background suppression, such as usage of the $\Sigma_{b}^{(*)}$
(or $B^{\ast 0}_{s2}$) tagging, more complex decay chains or the use of
control regions might be needed to perform the proposed searches at
significantly higher integrated luminosity. The presence of systematic
uncertainties might also risk reaching the whole region of theoretical
interest, with branching fractions $\gtrsim 10^{-5}$ for the channels of
interest and any value of $m_{\psi_{\rm{DS}}}$. Background modeling will of
course be particularly crucial when $m_{\psi_{\rm{DS}}}$ approaches the mass
of an existing SM resonance, since isolation requirements only reduce the
backgrounds from $b\bar{b}$ pairs, but do not grant that the remaining signal
candidates are not coming from exclusive backgrounds. More details about the
effect of systematic uncertainties are explained in A.

##  6 Conclusions

The B-mesogenesis mechanism predicts branching fractions of $b$-hadron decays
to dark baryons at relatively large values ($\sim 10^{-2}-10^{-6}$). The
closer the CP-violation sources are to the SM prediction, the higher those
branching fractions need to be in order to fulfill the baryon asymmetry of the
universe. In our study, we show that LHCb Upgrade will have the statistical
sensitivity to search for branching fractions at the $10^{-5}$ level or better
for the entire mass range of the dark baryon $\psi_{\rm{DS}}$. Multi-body
final states, such as
$B^{+}\rightarrow\psi_{\rm{DS}}(940)\Lambda_{c}(2595)^{+}$, show the best
performance at lower values of $m_{\psi_{\rm{DS}}}$ while two-track final
states are needed to reach higher masses, with $\Lambda_{b}^{0}$ decays
allowing to test the full allowed mass range. Systematic uncertainties could
play an important role in the sensitivity, and hence very precise background
modelling and suppression could be needed.

Together with a very precise measurement of the $B_{s}^{0}$ mixing parameters
and improved theoretical understanding of the model, the LHCb upgraded
experiment could arguably exclude completely, or confirm, B-mesogenesis as the
underlying mechanism for the baryon asymmetry of the universe as well as the
explanation for the Dark Matter problem.

###### Acknowledgements.

We would like to thank Miguel Escudero, Gonzalo Alonso-Álvarez, Gilly Elor,
Luca Silvestrini, Kristof De Bruyn, Lucia Grillo, and Alex Lenz for useful
discussions. This work has received financial support from Xunta de Galicia
(Centro singular de investigación de Galicia accreditation 2019-2022), by
European Union ERDF, and by the “María de Maeztu” Units of Excellence program
MDM-2016-0692 and the Spanish Research State Agency. In particular: the work
of X.C.V. is supported by MINECO (Spain) through the Ramón y Cajal program
RYC-2016-20073 and by XuntaGAL under the ED431F 2018/01 project.

##  Appendix A Effect of systematic uncertainties

The sensitivity results presented in this document make no assumptions on the
systematic uncertainties. However, as mentioned in Sec. 5, one would expect
them to limit the actual experimental analysis. We expect the dominant
systematic uncertainty to arise from the knowledge of the actual background
yield, which originates, among others, from the uncertainty in the $b\bar{b}$
cross section Aaij:2016avz . While any experimental analysis will do its best
to cure this uncertainty, this will eventually limit the precision to be
achieved. Potential strategies to mitigate this uncertainty include the use of
heavier $b$ resonances to tag the events of interest, as discussed in Sec.
3.5, the introduction of smarter selection criteria to discriminate the signal
against the background or the use of control regions that would provide a
handle to determine the actual background content in the signal region.

To provide a more quantitative estimate of the effect of the systematic
uncertainty in the background yield, we repeat the results of Sec. 4 for two
of our benchmark channels, with the same methodology but now under different
assumptions on the value of this uncertainty. These can found in Fig. 6 and
Fig. 7 for the $B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$ and
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ channels, respectively.
The main conclusion is that, in a pessimistic scenario, a systematic
uncertainty of $O(10)\%$ would make the search systematically limited already
with $\sim 1\,\rm{fb}^{-1}$, compromising the reach of all the region of
theoretical interest, with branching fractions $\gtrsim 10^{-5}$, for the
channels of interest and any value of $m_{\psi_{\rm{DS}}}$. This region would
be completely at reach with uncertainties of $O(1)\%$ or smaller. Finally,
bringing the uncertainty down to the level per-mile would increase the
systematic limit of the search to $\sim 15\,\rm{fb}^{-1}$.

Figure 6: Projected sensitivity in terms of branching fractions excluded at
95% CL for the $B^{0}\rightarrow\psi_{\rm{DS}}\Lambda(1520)$ channel. The
sensitivity is shown assuming integrated luminosities of 1, 15, 50 and 300
$\textrm{fb}^{-1}$. In each case, different assumptions are made regarding the
systematic uncertainty of the background yield, ranging from no uncertainty at
all to 10%. Figure 7: Same as Fig. 6, but now for the
$B^{+}\rightarrow\psi_{\rm{DS}}\Lambda_{c}(2595)^{+}$ channel.

## References

  * (1) G. Elor, M. Escudero, A. Nelson, Phys. Rev. D 99(3), 035031 (2019).  DOI 10.1103/PhysRevD.99.035031
  * (2) G. Alonso-Álvarez, G. Elor, M. Escudero, Phys. Rev. D 104(3), 035028 (2021).  DOI 10.1103/PhysRevD.104.035028
  * (3) Y.S. Amhis, et al., Eur. Phys. J. C 81(3), 226 (2021).  DOI 10.1140/epjc/s10052-020-8156-7.  arxiv:1909.12524. Updated results and plots available at https://hflav.web.cern.ch
  * (4) A. Lenz, G. Tetlalmatzi-Xolocotzi, JHEP 07, 177 (2020).  DOI 10.1007/JHEP07(2020)177
  * (5) H.M. Asatrian, H.H. Asatryan, A. Hovhannisyan, U. Nierste, S. Tumasyan, A. Yeghiazaryan, Phys. Rev. D 102(3), 033007 (2020).  DOI 10.1103/PhysRevD.102.033007
  * (6) R. Barate, et al., Eur. Phys. J. C 19, 213 (2001).  DOI 10.1007/s100520100612
  * (7) G. Elor, R. McGehee, Phys. Rev. D 103(3), 035005 (2021).  DOI 10.1103/PhysRevD.103.035005
  * (8) F. Elahi, G. Elor, R. McGehee (2021).  arxiv:2109.09751
  * (9) V. Chobanova, D. Martínez Santos, C. Prouve, M. Romero Lamas (2020).  arxiv:2012.02692
  * (10) M. Borsato, et al. (2021).  arxiv:2105.12668
  * (11) I. Bediaga, et al., Tech. Rep. CERN-LHCC-2013-021. LHCB-TDR-013 (2013).  URL https://cds.cern.ch/record/1624070
  * (12) I. Bediaga, et al., Tech. Rep. CERN-LHCC-2012-007 (2012).  URL https://cds.cern.ch/record/1443882
  * (13) R. Aaij, et al., Tech. rep. (2018).  URL https://cds.cern.ch/record/2636441.  ISBN 978-92-9083-494-6
  * (14) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C.O. Rasmussen, P.Z. Skands, Comput. Phys. Commun. 191, 159 (2015).  DOI 10.1016/j.cpc.2015.01.024
  * (15) R. Aaij, et al., Phys. Rev. Lett. 118(5), 052002 (2017).  DOI 10.1103/PhysRevLett.119.169901,10.1103/PhysRevLett.118.052002.  [Erratum: Phys. Rev. Lett.119,no.16,169901(2017)]
  * (16) R. Aaij, et al., Phys. Rev. D 104(3), 032005 (2021).  DOI 10.1103/PhysRevD.104.032005
  * (17) R. Aaij, et al., Tech. rep., CERN, Geneva (2017).  URL https://cds.cern.ch/record/2244311
  * (18) A.A. Alves Jr, et al., Tech. Rep. CERN-LHCC-2014-001. LHCB-TDR-015 (2014).  URL https://cds.cern.ch/record/1647400
  * (19) X. Cid Vidal, P. Ilten, J. Plews, B. Shuve, Y. Soreq, Phys. Rev. D 100(5), 053003 (2019).  DOI 10.1103/PhysRevD.100.053003
  * (20) X. Cid Vidal, Y. Tsai, J. Zurita, JHEP 01, 115 (2020).  DOI 10.1007/JHEP01(2020)115
  * (21) P.A. Zyla, et al., PTEP 2020(8), 083C01 (2020).  DOI 10.1093/ptep/ptaa104
  * (22) R. Aaij, et al., Phys. Lett. B 699, 330 (2011).  DOI 10.1016/j.physletb.2011.04.031
  * (23) D. Karlen, Computers in Physics 12(4), 380 (1998).  DOI 10.1063/1.168691.  URL https://aip.scitation.org/doi/abs/10.1063/1.168691
  * (24) I. Bediaga, et al., Tech. Rep. CERN-LHCC-2013-022. LHCB-TDR-014 (2013).  URL https://cds.cern.ch/record/1624074
  * (25) M. Adinolfi, et al., Eur. Phys. J. C 73, 2431 (2013).  DOI 10.1140/epjc/s10052-013-2431-9
  * (26) S. Stone, L. Zhang, Adv. High Energy Phys. 2014, 931257 (2014).  DOI 10.1155/2014/931257
  * (27) R. Aaij, et al., Phys. Rev. Lett. 122(1), 012001 (2019).  DOI 10.1103/PhysRevLett.122.012001
  * (28) R. Aaij, et al., JHEP 06, 129 (2020).  DOI 10.1007/JHEP06(2020)129
  * (29) T. Junk, Nucl. Instrum. Meth. A 434, 435 (1999).  DOI 10.1016/S0168-9002(99)00498-2
  * (30) A.L. Read, J. Phys. G28, 2693 (2002).  DOI 10.1088/0954-3899/28/10/313