# Time-reversal asymmetries in $\Lambda_{b}$ semileptonic decays

Chao-Qiang<EMAIL_ADDRESS>Xiang-Nan
<EMAIL_ADDRESS>and Chia-Wei<EMAIL_ADDRESS>School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute
for Advanced Study, UCAS, Hangzhou 310024, China
University of Chinese Academy of Sciences, 100190 Beijing, China

###### Abstract

We study the decays of $\Lambda_{b}\to\Lambda_{c}(\to
B_{n}f)\ell^{-}\overline{\nu}$ with $\ell=e,\mu,\tau$, where $B_{n}$ and $f$
are the daughter baryons and the rest of the particles in $\Lambda_{c}$
cascade decays, respectively. In particular, we examine the full angular
distributions with polarized $\Lambda_{b}$ and lepton mass effects, in which
the time-reversal asymmetries are identified. We concentrate on the decay
modes of $\Lambda_{b}\to\Lambda_{c}(\to pK^{-}\pi^{+})\ell^{-}\overline{\nu}$
to demonstrate their experimental feasibility. We show that the observables
associated with the time-reversal asymmetries are useful to search for new
physics as they vanish in the standard model. Explicitly, we find that they
are sensitive to the right-handed current from new physics, and possible to be
observed at LHCb.

##  I Introductions

The transitions of $b\to c\ell^{-}\overline{\nu}$ with $\ell=e,\mu,\tau$ have
raised great interest in both the theoretical and experimental aspects Bmeson
; BmesonEXP . In particular, the discrepancy of $R_{D^{(*)}}=\Gamma(\bar{B}\to
D^{(*)}\tau\overline{\nu})/\Gamma(\bar{B}\to D^{(*)}l\overline{\nu})$ with
$l=e,\mu$ has shown that the possible contributions from new physics (NP) can
be as large as ${\cal O}(10\%)$. Explicitly, we have that
$R_{D,D^{*}}^{\text{exp}}=(0.340\pm 0.030,0.295\pm 0.014)$ from the
experiments BmesonEXP ; pdg ; HFLAV:2019otj  and
$R_{D,D^{*}}^{\text{SM}}=(0.304\pm 0.003,0.259\pm 0.006)$ from the Lattice QCD
calculation LatticeBmeson , implying that NP can play a significant role. For
a review, one is referred to Ref. ReVonLUV .

On the other hand, LHCb has recently announced the baryonic version of the
ratio to be
$R_{\Lambda_{c}}=\Gamma(\Lambda_{b}\to\Lambda_{c}\tau^{-}\overline{\nu})/R_{\Lambda_{c}}=\Gamma(\Lambda_{b}\to\Lambda_{c}l^{-}\overline{\nu})=0.242\pm
0.026\pm 0.040\pm 0.059$ LbMeasurements tau , where the first and second
uncertainties are statistical and systematic, respectively, and the third one
comes from the normalization channel of
$\Lambda_{b}\to\Lambda_{c}\pi^{+}2\pi^{-}$. In contrast to $R_{D^{(*)}}$,
$R_{\Lambda_{c}}$ is found to be larger in theory, given as
$R_{\Lambda_{c}}=0.324\pm 0.004$ based on Lattice QCD Lattice . Such opposite
behavior indicates that there would be some theoretical errors, which have not
been properly considered. Thus, as a complementarity, it is useful to examine
the angular distributions Unpolarized ; Korner ; Tau Decay . In most of works
in the literature, $\Lambda_{b}$ is assumed to be unpolarized. However, it is
important to analyze the polarized cases, since the polarization fraction
$P_{b}$ is recently found to be around $3\%$ in proton-proton collisions at
centre-of-mass energies of 13 TeV polarizedEXP . We emphasize that with
$P_{b}\neq 0$, the time-reversal (TR) asymmetries can be observed without the
cascade decays of $\Lambda_{c}$ as we will show in this work. Moreover, the
value of $3\%$ is twice larger than ${\cal
B}(\Lambda_{c}\to\Lambda\pi^{+},pK_{S}^{0})$, and hence it is useful to study
the cases with $P_{b}\neq 0$ for probing the TR asymmetries.

The angular distribution of $\Lambda_{b}\to\Lambda_{c}(\to
pK_{S}^{0})\mu^{-}\overline{\nu}$ with polarized $\Lambda_{b}$ was first given
in Ref. polarized . In this work, we provide the full angular distributions of
$\Lambda_{b}\to\Lambda_{c}(\to B_{n}f)\ell^{-}\overline{\nu}$, where $B_{n}$
is the daughter baryon and $f$ stands for the rest of the daughter particles.
In contrast to those in the literature, we extend the study to the three-body
$\Lambda_{c}$ decays to include $\Lambda_{c}\to pK^{-}\pi^{+}$ and
$\Lambda_{c}\to\Lambda l^{+}\nu$. In particular, $\Lambda_{c}\to
pK^{-}\pi^{+}$ has a great advantage for the experimental detection, since all
the particles in the final states are charged.

In the standard model (SM), the TR asymmetries in
$\Lambda_{b}\to\Lambda_{c}(\to B_{n}f)\ell^{-}\overline{\nu}$ are zero due to
the absence of the weak phase in the $\Lambda_{b}\to\Lambda_{c}$ transition.
Clearly, a non-vanishing TR asymmetry indicates the existence of NP with a new
CP violating phase beyond the SM.

The layout of this work is given as follows. In Sec. II, we present the
angular distributions of the SM parameterized by the helicity amplitudes. In
Sec. III, we discuss the effects from NP, and show that they can be absorbed
by the redefinitions of the helicity amplitudes. In Sec. IV, we estimate the
TR asymmetries and their feasibility to be measured at LHCb. At last, we
conclude the study in Sec. V.

##  II Decay observables

In the SM, the amplitudes of $\Lambda_{b}\to\Lambda_{c}\ell^{-}\overline{\nu}$
are dominated by the weak interaction at tree level, given as

$\displaystyle\frac{G_{F}}{\sqrt{2}}V_{cb}g^{\mu\nu}\overline{u}_{\ell}{\gamma_{\mu}}(1-{\gamma_{5}})v\langle\Lambda_{c}|\overline{c}{\gamma_{\nu}}(1-{\gamma_{5}})b|\Lambda_{b}\rangle,$
(1)

where $G_{F}$ is the Fermi constant, $V_{cb}$ corresponds to the Cabibbo-
Kobayashi-Maskawa (CKM) matrix element, and $u_{\ell}$ and $v$ are the Dirac
spinors of charged leptons and antineutrinos, respectively. In this work, we
do not specify the flavors of (anti)neutrinos as they can not be distinguished
in experiments.

We further decompose the amplitudes by expanding the Minkowsiki metric,

$g^{\mu\nu}=\varepsilon_{t}^{\mu}(q)\varepsilon_{t}^{\nu}(q)-\sum_{\lambda=0,\pm}\varepsilon^{\mu}_{\lambda}(q)\varepsilon^{\ast\nu}_{\lambda}(q)\,,$
(2)

where $q=(q^{0},\vec{q}\,)$ and $\varepsilon$ are the 4-momentum and
polarization vector of the off-shell $W$ boson, respectively. The subscript in
$\varepsilon$ denotes the helicity, where $t$ indicates timelike while the
others spacelike. In particular, we have that

${\varepsilon}^{\mu}_{\pm}=\frac{1}{\sqrt{2}}(0,\pm
1,i,0)^{T}\,,\quad{\epsilon}_{0}^{\mu}=(0,0,0,-1)^{T}\,,\quad{\varepsilon}^{\mu}_{t}=(-1,0,0,0)^{T}\,,$
(3)

in the center of mass frame of $\ell^{-}\overline{\nu}$, which would be
referred to as the $\vec{q}$ frame in the following. Notice that the relative
phases between $\varepsilon$ are crucial as they interfere in the decay
distributions. In this work, they are fixed by the lowering operators, given
by

$\left(J_{x}-iJ_{y}\right)\epsilon_{1,0}=\sqrt{2}\epsilon_{0,-1}\,,$ (4)

where $J_{x,y}$ are the $SO(3)$ rotational generators. On the other hand, in
the center of mass frame of $\Lambda_{b}$ with $\vec{q}=-|\vec{q}\,|\hat{z}$,
which would be referred to as the $\Lambda_{b}$ frame, we have

${\epsilon}^{\mu}_{\pm}=\frac{1}{\sqrt{2}}(0,\mp
1,i,0)^{T}\,,\quad{\epsilon}_{0}^{\mu}=\frac{1}{\sqrt{q^{2}}}(-|\vec{q}\,|,0,0,q^{0})^{T}\,,\quad{\epsilon}^{\mu}_{t}=\frac{-1}{\sqrt{q^{2}}}q^{\mu},$
(5)

which are useful for the latter purpose.

Plugging Eq. (2) in Eq. (1), we have

$\frac{G_{F}}{\sqrt{2}}V_{cb}\left(L_{t}B_{t}-\sum_{\lambda=0,\pm}L_{\lambda}B_{\lambda}\right)\,,$
(6)

and

$B_{\lambda_{W}}=\epsilon_{\lambda_{W}}^{\ast\mu}\langle\Lambda_{c}|\overline{c}{\gamma_{\mu}}(1-{\gamma_{5}})b|\Lambda_{b}\rangle\,,\quad
L_{\lambda_{W}}=\epsilon_{\lambda_{W}}^{\mu}\overline{u}_{\ell}{\gamma_{\mu}}(1-{\gamma_{5}})v\,,$
(7)

with ${\lambda_{W}}=t,0$ and $\pm$. Note that $B_{\lambda_{W}}$ and
$L_{\lambda_{W}}$ depend on the polarizations of the baryons and leptons,
respectively. It is clear that in Eqs. (6) and (7), the amplitudes are
decomposed as the products of Lorentz scalars, describing
$\Lambda_{b}\to\Lambda_{c}W^{-*}$ ($B_{\lambda_{W}}$) and $W^{-*}\to
l^{-}\overline{\nu}$ ($L_{\lambda_{W}}$). A great advantage is that
$B_{\lambda_{W}}$ and $L_{\lambda_{W}}$ can be computed independently in the
$\Lambda_{b}$ and $\vec{q}$ frames, respectively, reducing the three-body
problems to the products of two-body ones.

To proceed further, we have to consider the polarizations of the baryons and
leptons. To this end, it is convenient to parameterized $B_{\lambda_{W}}$ as

$\displaystyle
B_{\lambda_{W}}=\epsilon^{\ast\mu}_{\lambda_{W}}\overline{u}_{c}$
$\displaystyle\left[\left(f_{1}(q^{2}){\gamma_{\mu}}-if_{2}(q^{2})\frac{\sigma_{\mu\nu}}{M_{b}}q^{\nu}+f_{3}(q^{2})\frac{q_{\mu}}{M_{b}}\right)\right.$
(8)
$\displaystyle~{}~{}~{}-\left.\left(g_{1}(q^{2}){\gamma_{\mu}}-ig_{2}(q^{2})\frac{\sigma_{\mu\nu}}{M_{b}}q^{\nu}+g_{3}(q^{2})\frac{q_{\mu}}{M_{b}}\right)\right]{\gamma_{5}}u_{b},$

where $f_{1,2,3}$ and $g_{1,2,3}$ represent the form factors, $M_{b,c}$ are
the masses of $\Lambda_{b,c}$, and
$\sigma_{\mu\nu}=i(\gamma_{\mu}\gamma_{\nu}-\gamma_{\nu}\gamma_{\mu})/2$. The
helicity amplitudes are calculated by

$\displaystyle
H_{\lambda_{c},\lambda_{W}}=B_{\lambda_{W}}\left(\lambda_{b}=\lambda_{c}-\lambda_{W},\lambda_{c},\vec{p}_{c}=-\vec{q}=|\vec{p}_{c}|\hat{z}\right)\,,$
(9)

where $\lambda_{b(c)}$ corresponds to the angular momentum (helicity) of
$\Lambda_{b(c)}$, $\vec{p}_{c}$ is the 3-momentum of $\Lambda_{c}$ in the
$\Lambda_{b}$ frame, and the conventions of the Dirac spinors are given in
Appendix A. Plugging Eq. (5) in Eq. (8), we obtain explicitly that

$\displaystyle H_{\pm\frac{1}{2}\pm
1}=\sqrt{2Q_{-}}\left(f_{1}+\frac{M_{+}}{M_{b}}f_{2}\right)\pm\sqrt{2Q_{+}}\left(-g_{1}+\frac{M_{-}}{M_{b}}g_{2}\right)\,,$
$\displaystyle
H_{\pm\frac{1}{2}0}=-\sqrt{\frac{Q_{-}}{q^{2}}}\left(M_{+}f_{1}+\frac{q^{2}}{M_{b}}f_{2}\right)\pm\sqrt{\frac{Q_{+}}{q^{2}}}\left(M_{-}g_{1}-\frac{q^{2}}{M_{b}}g_{2}\right)\,,$
$\displaystyle
H_{\pm\frac{1}{2}t}=-\sqrt{\frac{Q_{+}}{q^{2}}}\left(M_{-}f_{1}+\frac{q^{2}}{M_{b}}f_{3}\right)\pm\sqrt{\frac{Q_{-}}{q^{2}}}\left(M_{+}g_{1}-\frac{q^{2}}{M_{b}}g_{3}\right)\,,$
(10)

where $M_{\pm}=M_{b}\pm M_{c}$ and $Q_{\pm}=(M_{\pm})^{2}-q^{2}$. Note that
both the form factors and amplitudes depend on $q^{2}$.

On the other hand, the antineutrinos have positive helicities, and
$L_{\lambda_{W}}$ depends only on the helicity of $\ell^{-}$. By combining
Eqs. (3) and (7), we find

$h_{-}=-2\sqrt{2(q^{2}-m_{\ell}^{2})}\,,~{}~{}~{}h_{+}=\sqrt{\delta_{\ell}}h_{-}\,,~{}~{}~{}\delta_{\ell}=\frac{m_{\ell}^{2}}{2q^{2}}\,,$
(11)

where

$h_{\pm}=L_{-\frac{1}{2}\pm\frac{1}{2}}\left(\lambda_{e}=\pm\frac{1}{2},\vec{p}_{\ell}=-\vec{p}_{\nu}=|\vec{p}_{e}|\hat{z}\right)\,,$
(12)

with $L_{t}=L_{0}$ , and $\vec{p}_{\ell(\nu)}$ the 3-momentum of
$\ell^{-}(\overline{\nu})$ in the $\vec{q}$-frame .

The angular distributions of $\Lambda_{b}\to\Lambda_{c}(\to
B_{n}f)\ell^{-}\overline{\nu}$ can be obtained by piling up the Wigner-d
matrices of $d^{J}$, read as

$\displaystyle\frac{\partial^{6}\Gamma}{\partial
q^{2}\partial\cos\theta_{b}\partial\cos\theta_{c}\partial\cos\theta_{\ell}\partial\phi_{c}\partial\phi_{\ell}}\propto\frac{(q^{2}-m_{\ell}^{2})|\vec{p}_{c}|}{q^{2}}\sum_{\lambda_{\ell}\,,\lambda\,,\lambda_{b}}\rho_{\lambda_{b},\lambda_{b}}\left|A^{c}_{\lambda}h_{\lambda_{\ell}}\right|^{2}$
$\displaystyle\left|\sum_{\lambda_{c}\,,\lambda_{W}}(-1)^{J_{W}}H_{\lambda_{c}\,,\lambda_{W}}d^{\frac{1}{2}}(\theta_{b})^{\lambda_{b}}\,_{\lambda_{c}-\lambda_{W}}d^{\frac{1}{2}}(\theta_{c})^{\lambda_{c}}\,_{\lambda}d^{J_{W}}(\theta_{\ell})^{\lambda_{W}}\,_{\lambda_{\ell}-\frac{1}{2}}e^{i(\lambda_{c}\phi_{c}+\lambda_{\ell}\phi_{\ell})}\right|^{2},$
(13)

where $\rho_{\pm,\pm}=(1\pm P_{b})/2$, $\lambda_{(b,c,\ell)}=\pm 1/2$ ,
$|\vec{p}_{c}|=\sqrt{Q_{+}Q_{-}}/2M_{{\bf B}_{b}}$, the factor of
$(-1)^{J_{W}}$ comes from Eq. (6) along with $J_{W}=0~{}(1)$ for
$\lambda_{W}=t~{}(\pm,0)$, and $A_{\lambda}^{c}$ are associated with the up-
down asymmetries of $\Lambda_{c}\to B_{n}f$. Here, the definitions of the
angles can be found in FIG. 1, where $\theta_{b,c}$ and $\theta_{\ell}$ are
defined in the center of mass frames of $\Lambda_{b,c}$ and
$\ell^{-}\overline{\nu}$, respectively, while $\phi_{c,\ell}$ are the
azimuthal angles between the decay planes.

In the case of $\Lambda_{c}\to pK^{-}\pi^{+}~{}(\Lambda l^{+}\nu)$,
$A_{\lambda}^{c}$ depends on the 3-momentum of $p(\Lambda)$ and angles in
$K^{-}\pi^{+}$ $(l^{+}\nu)$ as well. However, we integrate out the dependence
for simplicity in this work. In addition, the cascade decays of $\tau^{-}$ can
be included by continually piling up the Wiger-d matrices inside Eq. (II). The
interested readers are referred to Ref. Tau Decay . Note that the overall
$q^{2}$ dependence in Eq. (II) can be cast in a more symmetric form by
recognizing $|\vec{p}_{\ell}|=(q^{2}-m_{\ell}^{2})/\sqrt{4q^{2}}$ in the
$\vec{q}$-frame.

Figure 1:  Definitions of the angles, where $B_{n}$ represents the daughter
baryon and $f$ the rest of the decay particles.

We expand the angular distributions as

$\displaystyle\begin{aligned} \frac{\partial^{6}\Gamma}{\partial
q^{2}\partial\vec{\Omega}}\propto\frac{(q^{2}-m_{\ell}^{2})^{2}|\vec{p}_{c}|}{q^{2}}\sum_{i=1}^{26}{\cal
X}_{i}(a_{\pm},b_{\pm},t_{\pm}){\cal P}_{i}(\alpha,P_{b}){\cal
D}_{i}(\vec{\Omega}),\end{aligned}$ (14)

where $\vec{\Omega}=(\theta_{b,c,\ell}\,,\phi_{c,\ell})$, and explicit forms
of ${\cal X}_{i}$, ${\cal P}_{i}$ and ${\cal D}_{i}$ can be found in Table 1,
where we have taken the abbreviations:

$a_{\pm}=H_{\pm\frac{1}{2},0}\,,~{}~{}~{}b_{\pm}=H_{\mp\frac{1}{2},\mp
1}\,,~{}~{}~{}t_{\pm}=H_{\pm\frac{1}{2},t}\,,~{}~{}~{}|\xi^{2}|=|\xi_{+}^{2}|+|\xi_{-}^{2}|\,,~{}~{}~{}|\xi_{\Delta}^{2}|=|\xi_{+}^{2}|-|\xi_{-}^{2}|\,,$
(15)

with $\xi=a,b,c$ and

$P_{2}=\frac{1}{2}\left(3\cos\theta_{\ell}-1\right)\,,~{}~{}~{}{\cal
X}_{i}{\cal D}_{i}=\text{Re}({\cal X}_{i}^{\phi})\text{Re}({\cal
D}_{i}^{\phi})-{\bf\text{Im}}({\cal X}_{i}^{\phi}){\bf\text{Im}}({\cal
D}_{i}^{\phi})\,.$ (16)

For instance, we have

$\displaystyle{\cal X}_{7}{\cal
D}_{7}=\frac{3}{\sqrt{2}}\left[\text{Re}\right.$
$\displaystyle\left(2\delta_{\ell}(t_{-}b_{+}^{*}-b_{-}t_{+}^{*})+(a_{-}b_{+}^{*}+b_{-}a_{+}^{*})\right)\cos\phi_{\ell}\sin\theta_{b}\sin\theta_{\ell}-$
(17)
$\displaystyle\left.{\bf\text{Im}}\left(2\delta_{\ell}(t_{-}b_{+}^{*}-b_{-}t_{+}^{*})+(a_{-}b_{+}^{*}+b_{-}a_{+}^{*})\right)\sin\phi_{\ell}\sin\theta_{b}\sin\theta_{\ell}\right]\,,$
$\displaystyle{\cal X}_{12}{\cal
D}_{12}=\frac{3}{\sqrt{2}}\left[\text{Re}\right.$
$\displaystyle\left(2\delta_{\ell}(b_{-}t_{-}^{*}-t_{+}b_{+}^{*})-(a_{+}b_{+}^{*}+b_{-}a_{-}^{*})\right)\cos\Phi\sin\theta_{c}\cos\theta_{c}-$
$\displaystyle\left.{\bf\text{Im}}\left(2\delta_{\ell}(b_{-}t_{-}^{*}-t_{+}b_{+}^{*})-(a_{+}b_{+}^{*}+b_{-}a_{-}^{*})\right)\sin\Phi\sin\theta_{c}\cos\theta_{c}\right]\,.$

Notice that $\xi_{\pm}$ are real in the SM, and any observations of nonzero
${\bf\text{Im}}({\cal X}_{i})$ would be a smoking gun of NP. For those ${\cal
D}_{i}^{\phi}$, which are independent of $\phi_{c,\ell}$, we simply have
${\cal X}_{i}{\cal D}_{i}={\cal X}_{i}^{\phi}{\cal D}_{i}^{\phi}$. The angular
distributions of
$\overline{\Lambda_{b}}\to\overline{\Lambda_{c}}(\to\overline{B_{n}}\bar{f})l^{+}\nu$
can be obtained directly by taking $\theta_{\ell}\to\pi-\theta_{\ell}$ and
$\alpha\to-\alpha$. In practice, $\delta_{l}$ can be taken as zero as an
excellent approximation in the SM, with which $t_{\pm}$ can be neglected as
well since they are always followed by $\delta_{\ell}$.

It is interesting to point out that, under the parity transformation, the
helicity amplitudes behave differently as

$a_{\pm}\to a_{\mp}\,,~{}~{}~{}b_{\pm}\to b_{\mp}\,,~{}~{}~{}t_{\pm}\to-
t_{\mp}\,,$ (18)

so that ${\cal X}_{7,12}^{R(I)}$ are parity even (odd). If $\Lambda_{b}$ is
unpolarized ($P_{b}=0$), it is clear that $\phi_{c}$ and $\phi_{\ell}$ can not
be measured separately. In this case, it is convenient to introduce a new set
of azimuthal coordinates as

$\displaystyle\Phi=\phi_{\ell}+\phi_{c}\,,\,\,\,\,\,\,\,\,\,$ $\displaystyle
0<\Phi<2\pi\,,$ (19)
$\displaystyle\Phi_{\Delta}=\frac{1}{2}(\phi_{\ell}-\phi_{c})\,,\,\,\,\,\,\,$
$\displaystyle-\pi<\Phi_{\Delta}<\pi\,.$

To obtain the unpolarized angular distributions from the polarized ones, one
can integrate over $\Phi_{\Delta}$ and $\cos\theta_{b}$, in which ${\cal
D}_{4\text{--}8,14\text{--}26}$ are zero. As a cross-check, we find that the
results are identical to those given in Ref. Korner .

Table 1: The angular distributions of $\Lambda_{b}\to\Lambda_{c}(\to
B_{n}f)\ell^{-}\overline{\nu}$ with $P_{2}=(3\cos^{2}\theta_{\ell}-1)/2$ and
the parameters of $a$, $b$ and $t$ defined in Eq. (15).

$i$ | ${\cal X}_{i}^{\phi}$ | $P_{i}$ | ${\cal D}_{i}^{\phi}$
---|---|---|---
1 | $(\delta_{\ell}+1)(|a|^{2}+|b|^{2})+3\delta_{\ell}|t|^{2}$ | $1$ | $1$
2 | $(2\delta_{\ell}-1)(|a|^{2}-\frac{1}{2}|b|^{2})$ | $1$ | $P_{2}$
3 | $-6\delta_{\ell}\left(\text{Re}\left(a_{+}t_{+}^{*}\right)+\text{Re}\left(a_{-}t_{-}^{*}\right)\right)+\frac{3}{2}|b_{\Delta}|^{2}$ | $1$ | $\cos\theta_{\ell}$
4 | $(\delta_{\ell}+1)(|a_{\Delta}|^{2}+|b_{\Delta}|^{2})+3\delta_{\ell}|t_{\Delta}|^{2}$ | $P_{b}$ | $\cos\theta_{b}$
5 | $(2\delta_{\ell}-1)(|a_{\Delta}|^{2}-\frac{1}{2}|b_{\Delta}|^{2})$ | $P_{b}$ | $\cos\theta_{b}P_{2}$
6 | $-6\delta_{\ell}\left[\text{Re}(a_{+}t_{+}^{*})-\text{Re}(a_{-}t_{-}^{*})\right]-\frac{3}{2}|b_{\Delta}|^{2}$ | $P_{b}$ | $\cos\theta_{b}\cos\theta_{\ell}$
7 | $\frac{3}{\sqrt{2}}\left[2\delta_{\ell}(t_{-}b_{+}^{*}-b_{-}t_{+}^{*})+(a_{-}b_{+}^{*}+b_{-}a_{+}^{*})\right]$ | $P_{b}$ | $e^{i\phi_{\ell}}\sin\theta_{b}\sin\theta_{\ell}$
8 | $\frac{3}{\sqrt{2}}(2\delta_{\ell}-1)(b_{-}a_{+}^{*}-a_{-}b_{+}^{*})$ | $P_{b}$ | $e^{i\phi_{\ell}}\sin\theta_{b}\sin\theta_{\ell}\cos\theta_{\ell}$
9 | $(\delta_{\ell}+1)(|a_{\Delta}|^{2}-|b_{\Delta}|^{2})+3\delta_{\ell}|t_{\Delta}|^{2}$ | $\alpha$ | $\cos\theta_{c}$
10 | $(2\delta_{\ell}-1)(|a_{\Delta}|^{2}+\frac{1}{2}|b_{\Delta}|^{2})$ | $\alpha$ | $\cos\theta_{c}P_{2}$
11 | $6\delta_{\ell}\left(\text{Re}\left(a_{-}t_{-}^{*}\right)-\text{Re}\left(a_{+}t_{+}^{*}\right)\right)-\frac{3}{2}|b|^{2}$ | $\alpha$ | $\cos\theta_{c}\cos\theta_{\ell}$
12 | $\frac{3}{\sqrt{2}}\left[2\delta_{\ell}(b_{-}t_{-}^{*}-t_{+}b_{+}^{*})-(a_{+}b_{+}^{*}+b_{-}a_{-}^{*})\right]$ | $\alpha$ | $e^{i\Phi}\sin\theta_{c}\cos\theta_{c}$
13 | $\frac{3}{\sqrt{2}}\left[2\delta_{\ell}(a_{+}b_{+}^{*}-b_{-}a_{-}^{*})+(b_{-}a_{-}^{*}-a_{+}b_{+}^{*})\right]$ | $\alpha$ | $e^{i\Phi}\sin\theta_{c}\cos\theta_{c}\cos\theta_{\ell}$
14 | $(2\delta_{\ell}+1)a_{+}a_{-}^{*}$ | $P_{b}\alpha$ | $\sin\theta_{b}\sin\theta_{c}P_{2}$
15 | $2\delta_{\ell}(|a|^{2}+|b|^{2})-(|a|^{2}+1/2|b|^{2})$ | $P_{b}\alpha$ | $\cos\theta_{b}\cos\theta_{c}P_{2}$
16 | $\delta_{\ell}(|a|^{2}-|b|^{2}+3|t|^{2})+(|a|^{2}-|b|^{2})$ | $P_{b}\alpha$ | $\cos\theta_{b}\cos\theta_{c}$
17 | $-6\delta_{\ell}\left(\text{Re}(a_{+}t_{+}^{*})+\text{Re}(a_{-}t_{-}^{*})\right)-\frac{3}{2}|b_{\Delta}|^{2}$ | $P_{b}\alpha$ | $\cos\theta_{b}\cos\theta_{c}\cos\theta_{\ell}$
18 | $(2-4\delta_{\ell})a_{-}a_{+}^{*}$ | $P_{b}\alpha$ | $e^{-i\phi_{c}}\sin\theta_{b}\sin\theta_{c}P_{2}$
19 | $(1-2\delta_{\ell})b_{-}b_{+}^{*}$ | $P_{b}\alpha$ | $e^{i(\phi_{c}+2\phi_{\ell})}\sin\theta_{b}\sin\theta_{c}P_{2}$
20 | $-2\delta_{\ell}(a_{+}a_{-}^{*}+3t_{+}t_{-}^{*})-2a_{+}a_{-}^{*}$ | $P_{b}\alpha$ | $e^{i\phi_{c}}\sin\theta_{b}\sin\theta_{c}$
21 | $(2\delta_{\ell}-1)b_{-}b_{+}^{*}$ | $P_{b}\alpha$ | $e^{i(\phi_{c}+2\phi_{\ell})}\sin\theta_{b}\sin\theta_{c}$
22 | $6\delta_{\ell}(a_{+}t_{-}^{*}+t_{+}a_{-}^{*})$ | $P_{b}\alpha$ | $e^{i\phi_{c}}\sin\theta_{b}\sin\theta_{c}\cos\theta_{\ell}$
23 | $\frac{3}{\sqrt{2}}\left[(b_{-}a_{+}^{*}-a_{-}b_{+}^{*})-2\delta_{\ell}(b_{-}t_{+}^{*}+t_{-}b_{+}^{*})\right]$ | $P_{b}\alpha$ | $e^{i\phi_{\ell}}\sin\theta_{b}\sin\theta_{\ell}\cos\theta_{c}$
24 | $\frac{3}{\sqrt{2}}\left[(b_{-}a_{-}^{*}-a_{+}b_{+}^{*})-2\delta_{\ell}(b_{-}t_{-}^{*}+t_{+}b_{+}^{*})\right]$ | $P_{b}\alpha$ | $e^{i\Phi}\sin\theta_{c}\sin\theta_{\ell}\cos\theta_{b}$
25 | $\frac{3}{\sqrt{2}}\left[2\delta_{\ell}(a_{-}b_{+}^{*}+b_{-}a_{+}^{*})-(a_{-}b_{+}^{*}+b_{-}a_{+}^{*})\right]$ | $P_{b}\alpha$ | $e^{i\phi_{\ell}}\sin\theta_{b}\sin\theta_{\ell}\cos\theta_{c}\cos\theta_{\ell}$
26 | $\frac{3}{\sqrt{2}}\left[2\delta_{\ell}(b_{-}a_{-}^{*}-a_{+}b_{+}^{*})-(a_{+}b_{+}^{*}+b_{-}a_{-}^{*})\right]$ | $P_{b}\alpha$ | $e^{i\Phi}\sin\theta_{c}\sin\theta_{\ell}\cos\theta_{b}\cos\theta_{\ell}$

With Table 1, one can construct several observables in a model independent
way. The simplest ones would be the partial and total decay widths, read as

$\displaystyle\frac{d\Gamma}{dq^{2}}=\zeta(q^{2}){\cal
X}_{1}\,,~{}~{}~{}~{}\zeta(q^{2})=\frac{G_{F}^{2}}{24\pi^{3}}V_{cb}^{2}\frac{(q^{2}-M_{\ell}^{2})^{2}|\vec{p}_{c}|}{8M_{{\bf
B}_{b}}^{2}q^{2}}\,,$ (20)

and

$\displaystyle\Gamma=\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta(q^{2}){\cal
X}_{1}dq^{2}\,,$ (21)

respectively. It shall be clear that
$\Gamma(\Lambda_{b}\to\Lambda_{c}\ell^{-}\overline{\nu})$ is independent of
$\Lambda_{c}\to B_{n}f$. Likewise, there are several angular observables that
can be defined independent of $\Lambda_{c}\to B_{n}f$, and it is reasonable to
measure them separately as they do not suffer from the smallness of ${\cal
B}(\Lambda_{c}\to B_{n}f)$. In fact, the angular distributions without cascade
decays can be obtained straightforwardly by integrating over
$(\cos\theta_{c},\phi_{c})$, resulting in

$\displaystyle\begin{aligned} \frac{\partial^{4}\Gamma}{\partial
q^{2}\partial\cos\theta_{b}\cos\theta_{\ell}\partial\phi_{\ell}}\propto\frac{(q^{2}-m_{\ell}^{2})^{2}|\vec{p}_{c}|}{q^{2}}\sum_{i=1}^{8}{\cal
X}_{i}(a_{\pm},b_{\pm},t_{\pm}){\cal P}_{i}(0,P_{b}){\cal
D}_{i}(\vec{\Omega}),\end{aligned}$ (22)

which is clearly independent of $\alpha$. As a cross-check, we find that Eq.
(22) reduces to the ones given by Ref. KornerNeutron  with an appropriate
substitution.

There are a couple of quantities that deserve a closer look. The forward-
backward asymmetries for $W^{-*}\to l^{-}\overline{\nu}$ and $\Lambda_{c}\to
B_{n}f$ are defined as

$\displaystyle
A_{FB}=2\left(\int^{1}_{0}-\int^{0}_{-1}\right)\Gamma_{\theta_{\ell}}d\cos\theta_{\ell}=\frac{1}{\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta(q^{2}){\cal
X}_{3}dq^{2}\,,$ $\displaystyle
A_{PL}=\frac{2}{\alpha}\left(\int^{1}_{0}-\int^{0}_{-1}\right)\Gamma_{\theta_{c}}d\cos\theta_{c}=\frac{1}{\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta(q^{2}){\cal
X}_{9}dq^{2}\,,$ (23)

where we have adopted the shorthand notation,

$\Gamma_{\Omega}=\frac{1}{\Gamma}\frac{\partial\Gamma}{\partial\vec{\Omega}}\,.$
(24)

The up-down asymmetries $A_{UD}$, on the other hand, are given by

$A_{UD}=\frac{2}{P_{b}}\left(\int^{1}_{0}-\int^{0}_{-1}\right)\Gamma_{\theta_{b}}d\cos\theta_{b}=\frac{1}{\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta(q^{2}){\cal
X}_{4}dq^{2}\,,$ (25)

which require $P_{b}\neq 0$ for an experimental measurement.

Here, it is an appropriate place to revisit $A_{\lambda}^{c}$ in Eq. (II)
explicitly. We have

$\left|A_{\pm\frac{1}{2}}^{c}\right|^{2}=\frac{1}{2}\left(1\pm\alpha\right)\,,$
(26)

where $\alpha$ are the up-down asymmetries of $\Lambda_{c}\to B_{n}f$, with
the experimental values given by pdg

$\displaystyle\begin{aligned}
\alpha(\Lambda_{c}\to\Lambda\pi^{+},&\Sigma^{0}\pi^{+},\Sigma^{+}\pi^{0},pK^{0}_{s})=\\\
&(-0.84\pm 0.09,-0.55\pm 0.11,-0.73\pm 0.18,0.2\pm 0.5),\end{aligned}$ (27)

Similarly, for the three-body $\Lambda_{c}$ decays, we have
$\alpha=A_{UD}^{c}$, where $A_{UD}^{c}$ are defined by substituting
$\Lambda_{c}\to B_{n}f$ for $\Lambda_{b}\to\Lambda_{c}\ell^{-}\overline{\nu}$
in Eq. (25). In particular, $\alpha$ are found to be

$\alpha(\Lambda_{c}\to pK^{-}\pi^{+},\Lambda l^{+}\nu)=(0.89\pm 0.10,-0.32),$
(28)

from the $SU(3)_{F}$ analysis UP-down  and light-front quark model lightfront
, respectively.

The azimuthal angles are closely related to the triple product asymmetries,
which flip signs under TR transformation TimeReversal . To probe them, we
define

$\displaystyle{\cal
T}_{\ell}=\frac{1}{P_{b}}\left(\int^{\pi}_{0}-\int^{\pi}_{2\pi}\right)\Gamma_{\phi_{\ell}}d\phi_{\ell}=-\frac{\pi^{2}}{8\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta{\bf\text{Im}}({\cal
X}_{7}^{\phi})dq^{2}\,,$ (29) $\displaystyle{\cal
T}_{c}=\frac{1}{\alpha}\left[\left(\int^{\pi}_{0}-\int^{\pi}_{2\pi}\right)d\Phi\right]\left[\left(\int_{0}^{1}-\int^{0}_{-1}\right)d\cos\theta_{c}\right]\Gamma_{\phi_{\ell},\theta_{c}}=-\frac{2}{3\pi\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta{\bf\text{Im}}({\cal
X}_{12}^{\phi})dq^{2}\,,$

which are proportional to the complex phases of $\xi_{\pm}$, and vanish
without NP. Comparing to the direct CP asymmetries, TR asymmetries do not
require strong phases, which are great advantages to probe CP violation as
strong phases are absent in the semileptonic decays. Note that one can also
construct other TR asymmetries from ${\cal X}_{8,13,18\text{--}26}$.

##  III Contributions from possible new physics

Let us consider the dimension-6 effective Hamiltonian from NP with left-handed
neutrinos, read as where $P_{R,L}=(1\pm\gamma_{5})$, $C_{S,P,R,L,T}$ are the
Wilson coefficients, which are complex and depend on the lepton flavors in
general, and $N$ in the superscript indicates that NP is considered. The
effects of $C_{S,P,R,L}$ can be absorbed by redefining the amplitudes as

$\displaystyle
a_{\pm}^{N}=\left(1+C_{L}\right)a_{\pm}+C_{R}a_{\mp}\,,~{}~{}~{}b^{N}_{\pm}=\left(1+C_{L}\right)b_{\pm}+C_{R}b_{\mp}\,,$
$\displaystyle
t_{\pm}^{N}=\left(1+C_{L}\right)t_{\pm}+C_{R}t_{\mp}-\frac{\sqrt{Q_{+}q^{2}}}{m_{\ell}}C_{S}f_{s}\pm\frac{\sqrt{Q_{-}q^{2}}}{m_{\ell}}C_{P}g_{p}\,,$
(30)

where $f_{s}$ and $g_{p}$ are defined by

$\langle\Lambda_{c}|\overline{c}\left(1+\gamma_{5}\right)b|\Lambda_{b}\rangle=\overline{u}_{c}\left(f_{s}+g_{p}\gamma_{5}\right)u_{b}\,.$
(31)

The derivations can be found in Appendix B. Note that in Eq. (III),
$\xi_{\pm}$ are calculated within the SM given in Eq. (II). The angular
distributions can be easily obtained by substituting $\xi_{\pm}^{N}$ for
$\xi_{\pm}$ in Table 1 . In the case of $C_{R,S,P}=0$, the effect of $C_{L}$
can be absorbed by redefining $V_{cb}$ as $V_{cb}(1+C_{L})$, leaving the
angular distributions unaltered. Therefore, in the following, we would simply
take $C_{L}=0$.

Let us first consider the case that $C_{R}\neq 0$ with $C_{S,P}=0$. For the
total decay widths, $C_{R}$ would be polluted by the uncertainties of the form
factors. However, we can utilize that $\xi_{\pm}$ are real, whereas $C_{R}$
can be complex in general. Plugging Eq. (III) in Eq. (17), we arrive

$\displaystyle{\bf\text{Im}}\left({\cal X}_{7}^{\phi
N}\right)=3\sqrt{2}\,{\bf\text{Im}}(C_{R})\left(a_{+}b_{+}-a_{-}b_{-}\right)\,,$
$\displaystyle{\bf\text{Im}}\left({\cal X}_{12}^{\phi
N}\right)=3\sqrt{2}\,{\bf\text{Im}}(C_{R})\left(a_{+}b_{-}-a_{-}b_{+}\right)\,,$
(32)

where we have taken $\xi_{\pm}$ as real, calculated by Eq. (II).

On the other hand, the effects of $C_{S,P}$ are largely enhanced by the
smallness of the lepton quark masses when $q^{2}/m_{\ell}^{2}\gg 1$ .
Therefore, measuring $t_{\pm}$ in high $q^{2}$ regions would be useful to
constrain the values of $C_{S,P}$. To diminish the uncertainties from the form
factors, one can examine the complex phases, given by

$\displaystyle{\bf\text{Im}}\left({\cal X}_{7}^{\phi
N}\right)=\frac{3}{\sqrt{2}}m_{\ell}\left[-{\bf\text{Im}}(C_{S})\sqrt{\frac{Q_{+}}{q^{2}}}f_{s}\left(b_{+}+b_{-}\right)+{\bf\text{Im}}(C_{P})\sqrt{\frac{Q_{-}}{q^{2}}}g_{p}\left(b_{-}-b_{+}\right)\right]\,,$
$\displaystyle{\bf\text{Im}}\left({\cal X}_{12}^{\phi
N}\right)=\frac{3}{\sqrt{2}}m_{\ell}\left[{\bf\text{Im}}(C_{S})\sqrt{\frac{Q_{+}}{q^{2}}}f_{s}\left(b_{+}+b_{-}\right)+{\bf\text{Im}}(C_{P})\sqrt{\frac{Q_{-}}{q^{2}}}g_{p}\left(b_{-}-b_{+}\right)\right]\,,$
(33)

where we have taken $C_{R}=0$. By collecting Eqs. (III) and (III), the net
effects of NP on ${\cal T}_{\ell,c}$ are summarized as follows

$\displaystyle{\cal T}_{\ell}$ $\displaystyle=$
$\displaystyle-\frac{3\pi^{2}}{8\sqrt{2}}\left({\bf\text{Im}}(C_{R}){\cal
Y}_{R}-{\bf\text{Im}}(C_{S}){\cal Y}_{S}+{\bf\text{Im}}(C_{P}){\cal
Y}_{P}\right)\,,$ $\displaystyle{\cal T}_{c}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{2}}{\pi}\left({\bf\text{Im}}(C_{R}){\cal
Y}^{\prime}_{R}+{\bf\text{Im}}(C_{S}){\cal Y}_{S}+{\bf\text{Im}}(C_{P}){\cal
Y}_{P}\right)\,,$ (34)

where

$\displaystyle{\cal Y}_{R}^{(\prime)}$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}2\zeta\left(a_{+}b_{+(-)}-a_{-}b_{-(+)}\right)dq^{2}\,,$
$\displaystyle{\cal Y}_{S}$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta
m_{\ell}\sqrt{\frac{Q_{+}}{q^{2}}}f_{s}\left(b_{+}+b_{-}\right)dq^{2}\,,$
$\displaystyle{\cal Y}_{P}$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma}\int^{(M_{b}-M_{c})^{2}}_{m_{\ell}^{2}}\zeta
m_{\ell}\sqrt{\frac{Q_{+}}{q^{2}}}g_{p}\left(b_{-}-b_{+}\right)dq^{2}\,.$ (35)

Notice that $\Gamma$ also depends on $C_{R,S,P}$. However, in this work, we
take $C_{R,S,P}$ as zero in $\Gamma$ as a first order approximation, and
therefore ${\cal Y}^{(\prime)}_{R,S,P}$ can be computed once the form factors
are given.

To examine TR asymmetries in the experiments, we define

$\displaystyle\Delta N_{\ell}\equiv$ $\displaystyle
N(\pi>\phi_{\ell}>0)-N(2\pi>\phi_{\ell}>\pi)=\epsilon N_{\Lambda_{b}}{\cal
T}_{\ell}{\cal B}(\Lambda_{b}\to\Lambda_{c}^{+}\ell^{-}\nu)P_{b}\,,$
$\displaystyle\Delta N_{f}\equiv$ $\displaystyle
N(\pi>\phi_{c}>0,\cos\theta_{c}>0)+N(2\pi>\phi_{c}>\pi,\cos\theta_{c}<0)$
$\displaystyle~{}~{}~{}~{}-N(\pi>\phi_{c}>0,\cos\theta_{c}<0)-N(2\pi>\phi_{c}>\pi,\cos\theta_{c}>0)$
$\displaystyle=$ $\displaystyle\epsilon N_{\Lambda_{b}}{\cal T}_{c}{\cal
B}(\Lambda_{b}\to\Lambda_{c}^{+}\ell^{-}\nu){\cal B}(\Lambda_{c}\to
B_{n}f)\alpha(\Lambda_{c}\to B_{n}f)\,,$ (36)

which hold at ${N\to\infty}$ with $N$ the number of the observed events, where
$N_{\Lambda_{b}}$ is the numbers of $\Lambda_{b}$ in experiments, and
$\epsilon$ is the efficiency for an experimental reconstruction. To reduce the
statistical uncertainties in $\Delta N_{f}$, we can sum over the decay modes,
given as

$\Delta N_{c}=\sum_{f}\left|N_{f}\right|=\epsilon N_{\Lambda_{b}}{\cal
T}_{c}{\cal B}(\Lambda_{b}\to\Lambda_{c}^{+}\ell^{-}\nu)\sum_{f}\left|{\cal
B}(\Lambda_{c}\to B_{n}f)\alpha(\Lambda_{c}\to B_{n}f)\right|\,.$ (37)

As $\Delta N_{{\ell,c}}$ are proportional to ${\cal T}_{\ell,c}$, a nonzero
value of $\Delta N_{{\ell}}$ or $\Delta N_{{c}}$ would be a smoking gun of NP.

The full angular distributions including the tensor operator are given in
Appendix C. For simplicity, we take $C_{T}=0$ in the numerical analysis, as
they can not be reduced to the form of Eq. (6), which breaks the angular
analysis. In addition, the tensor operator is closely related to the scalar
ones by the Fierz transformation in the leptoquak effective field theory
Leptoquark .

##  IV Numerical results

As mentioned in Sec. III, nonzero signals of ${\cal T}_{\ell,c}$ can be clear
evidences of CP violation from NP with $t^{N}_{\pm}$ largely enhanced in high
$q^{2}$ regions by $C_{S,P}$. In Eq. (III), ${\bf\text{Im}}(C_{R,S,P})$ are in
general free parameters of NP, whereas ${\cal Y}_{R,S,P}^{(\prime)}$ can be
computed by the form factors. In this work, we utilize the homogeneous bag
model to estimate the form factors Liu:2022pdk , which agree well with those
by the Lattice QCD calculation as well as the heavy quark symmetry Lattice .

The values of ${\cal Y}_{R,S,P}^{(\prime)}$ are listed in Table 2, from which
one can see that ${\cal Y}_{R}$ and ${\cal Y}_{R}^{\prime}$ are both sizable
and free of flavors, giving us a good opportunity to examine Im$(C_{R})$. In
contrast, the values of ${\cal Y}_{S,P}$ for $\ell=e$ and $\mu$ are suppressed
due to $m_{\ell}$. For $\ell=\tau$, they are still 3 times smaller than ${\cal
Y}_{R}^{(\prime)}$. Hence, Im$(C_{S,P})$ are much more difficult to be
observed in the experiments comparing to Im$(C_{R})$.

Table 2: Parameters defined in Eq. (III), where the uncertainties come from the model calculation.  $\ell$ | ${\cal Y}_{R}$ | ${\cal Y}_{R}^{\prime}$ | ${\cal Y}_{S}$ | ${\cal Y}_{P}$
---|---|---|---|---
$e$ | $-0.238(11)$ | $0.703(13)$ | $<10^{-4}$ | $<10^{-4}$
$\mu$ | $-0.237(12)$ | $0.701(13)$ | $-0.0037(2)$ | $0.0136(2)$
$\tau$ | $-0.149(6)$ | $0.438(8)$ | $-0.039(2)$ | $0.180(2)$

To explain the excesses of $R_{D^{(*)}}$, $C_{R,S,P}$ are found to be tiny for
$\ell=l$, but fortunately, $C_{R}=\pm 0.42(7)i$ is huge for $\ell=\tau$
Leptoquark . We have plotted $\zeta(\partial{\cal Y}_{R}^{(\prime)}/\partial
q^{2})$ for $\ell=\tau$ in FIG. 2, where the bands represent the uncertainties
from the form factors. One can see that the idea $q^{2}$ region to probe the
asymmetries lies around $q^{2}=9$ GeV2. Finally, putting the values of ${\cal
Y}_{R}^{(\prime)}$ and $C_{R}=\pm 0.42(7)i$ in Eq. (III), we find that

${\cal T}_{\ell}=\pm 0.16(3)\,,~{}~{}~{}~{}{\cal T}_{c}=\pm 0.08(2)\,,$ (42)

for $\Lambda_{b}\to\Lambda_{c}\tau^{-}\overline{\nu}$. Notice that the signs
are irrelevant for searching evidences of NP as long as they are nonzero. To
estimate the results at LHCb Run2, we take $N_{\Lambda_{b}}=5\times 10^{9}$,
$P_{b}=0.03$, $\varepsilon=10^{-4}$, and

$\sum_{f}\left|{\cal B}(\Lambda_{c}\to B_{n}f)\alpha(\Lambda_{c}\to
B_{n}f)\right|=6\times 10^{-2}\,,$ (43)

resulting in that $|\Delta N_{\ell}|\approx 50$ and $\Delta N_{c}\approx 20$
for $\ell=\tau$, which are large and ready to be measured. Here, Eq. (43) is
derived by crunching up the numbers in Eqs. (27) and (28).

Figure 2:  $q^{2}$ dependence of $\zeta(\partial{\cal
Y}_{R}^{(\prime)}/\partial q^{2})$ in
$\Lambda_{b}\to\Lambda_{c}\tau^{-}\overline{\nu}$, where the bands represent
the uncertainties caused by the form factors.

To probe the effects of the scalar operators, we find that ${\cal X}_{6}$,
which can be understood as a combination of $A_{FB}$ and $A_{UD}$, is
sensitive to $C_{S,P}$. The results are plotted in FIG. 3, where we have taken
$C_{P}=C_{S}$. In the region of 9 GeV${}^{2}<q^{2}<$10 GeV2, ${\cal X}_{6}$
can be enhanced largely. In particular, it is twice larger with $C_{S}=0.2$ in
comparison to that in the SM.

Figure 3:  $q^{2}$-dependence of ${\cal X}_{6}\zeta/\Gamma$ with $\ell=\tau$
for $C_{S}=0,0.1$ and $0.2$, respectively.

##  V Summary

Based on the helicity formalism, we have given the full angular distributions
of $\Lambda_{b}\to\Lambda_{c}(\to B_{n}f)\ell^{-}\overline{\nu}$. In
particular, we have identified TR violating terms, which vanish in the SM due
to the lack of relative complex phases. Since strong phases are not required
in these TR violating observables in contrast to the direct CP asymmetries,
they can be reliably calculated. The angular distributions have been given
explicitly with the helicity amplitudes in Table 1. We have cross-checked our
results with those in Refs. Korner ; KornerNeutron , and found that they are
consistent. Note that our results can be easily applied to
$\Xi_{b}\to\Xi_{c}(\to B_{n}f)\ell^{-}\overline{\nu}$ with trivial
modifications.

Notably, the effects of NP can be absorbed by redefining the helicity
amplitudes as demonstrated in Eq. (III) with $\xi_{\pm}$ calculated in the SM.
We recommend the experiments to measure the TR violating observables of ${\cal
T}_{\ell,c}$ defined in Eq. (29) for searching NP as they vanish in the SM. To
compare with the experiments, $\Delta N_{\ell,c}$ have been defined by the
numbers of the observe events, which are proportional to ${\cal T}_{\ell,c}$.
Based on $C_{R}=0.42(7)i$ for $\ell=\tau$, we have obtained that $|\Delta
N_{\ell,c}|\approx 50,20$ at LHCb Run2, which are sufficient for measurements.
On the other hand, we have pointed out that ${\cal X}_{6}$ is sensitive to
$C_{S,P}$ for $\Lambda_{b}\to\Lambda_{c}\tau^{-}\overline{\nu}$, which can be
largely enhanced in the high $q^{2}$ region.

###### Acknowledgements.

We would like to thank Jiabao Zhang for valuable discussions.

##  Appendix A Dirac spinors

In this work, we choose the fermion and antifermion in $p\hat{z}$ and
$-p\hat{z}$ directions, respectively. We have that

$u_{+}=\left(\begin{array}[]{c}E_{+}\\\ 0\\\ E_{-}\\\
0\end{array}\right)\,,\quad u_{-}=\left(\begin{array}[]{c}0\\\ E_{+}\\\ 0\\\
-E_{-}\end{array}\right)\,,\quad v_{+}=\left(\begin{array}[]{c}E_{-}\\\ 0\\\
-E_{+}\\\ 0\end{array}\right)\,,\quad v_{-}=\left(\begin{array}[]{c}0\\\
-E_{-}\\\ 0\\\ -E_{+}\end{array}\right)\,,$ (60)

with $``\pm"$ denoting the helicities, $E_{\pm}=\sqrt{E\pm m}$, and $m$ the
particle mass. Notice that the relative signs are crucial, fixed by the
relations

$u_{-}=L^{u}_{z}R_{y}(\pi)(L_{z}^{u})^{-1}u_{+}\,,~{}~{}~{}v_{-}=(L_{z}^{v})^{-1}R_{y}(\pi)L_{z}^{v}v_{+}\,,$
(61)

where $R_{y}$ is the rotation matrix toward $\hat{y}$, and $L_{z}^{u,v}$ are
the Lorentz boost operators toward $\hat{z}$. Here, $L_{z}^{u,v}$ are taken in
a way such that $L_{z}^{u-1}u_{\pm}$ and $L_{z}^{v}v_{\pm}$ are at rest.

##  Appendix B Contributions from scalar and current operators

We now scratch the derivations of Eq. (III). The contributions of $C_{L}$ are
given by

$\xi_{\pm}\to(1+C_{L})\xi_{\pm}\,,$ (62)

due to the same coupling in the SM. Furthermore, those of $C_{R}$ can be
obtained straightforwardly by

$\xi_{\pm}\to\xi_{\pm}+C_{R}\xi_{\mp}\,,$ (63)

as $\overline{c}\gamma^{\mu}P_{L,R}b$ are related by the parity.

The scalar operators contribute to the amplitudes as

$\displaystyle\frac{G_{f}}{\sqrt{2}}V_{cb}L_{t}^{N}B_{t}^{N}\,,~{}~{}~{}L_{t}^{N}={\cal
C}\overline{u}_{\ell}P_{L}v_{\nu}\,,~{}~{}~{}B_{t}^{N}={\cal
C}^{-1}\overline{u}_{c}\left(C_{S}f_{s}+C_{P}g_{p}\gamma_{5}\right)u_{b}\,,$
(64)

where ${\cal C}$ is a constant, Clearly, $B_{t}^{N}$ and $L_{t}^{N}$ can be
viewed as the transitions of $\Lambda_{b}\to\Lambda_{c}P^{-}$ and $P^{-}\to
l^{-}\overline{\nu}$, respectively, with $P^{-}$ an effective particle from
NP. As $P^{-}$ is spinless, $L_{t}^{N}B_{t}^{N}$ is only related to
$L_{t}B_{t}$ in Eq. (6). By adjusting ${\cal C}$ such that $L_{t}^{N}=L_{t}$,
we arrive

$t_{\pm}^{N}=t_{\pm}-\frac{\sqrt{Q_{+}q^{2}}}{m_{\ell}}C_{S}f_{s}\pm\frac{\sqrt{Q_{-}q^{2}}}{m_{\ell}}C_{P}g_{p}\,.$
(65)

By collecting Eqs. (62), (63) and (65), we obtain Eq. (III). It is interesting
to see that all the contributions can be encapsulated in $\xi_{\pm}$, which
already exist in the SM.

##  Appendix C Contributions from tensor operator

To cooperate the tensor operator with the helicity amplitudes, by utilizing
Eq. (2) we decompose the products of the Minkowsiki metric as

$\displaystyle\left(g^{\mu\nu}g^{\mu^{\prime}\nu^{\prime}}-g^{\mu^{\prime}\nu}g^{\mu\nu^{\prime}}\right)=\sum_{\lambda}\left(-V_{1}^{\mu\mu^{\prime}}(\lambda)V_{1}^{\nu\nu^{\prime}*}(\lambda)+V_{2}^{\mu\mu^{\prime}}(\lambda)V_{2}^{\nu\nu^{\prime}*}(\lambda)\right)\,,$
(66)

with

$\displaystyle
V_{1}^{\mu\mu^{\prime}}(\lambda)=\left(\varepsilon_{t}^{\mu}\varepsilon^{\mu^{\prime}}_{\lambda}-\varepsilon_{t}^{\mu^{\prime}}\varepsilon^{\mu}_{\lambda}\right)\,,~{}~{}~{}V_{2}^{\mu\mu^{\prime}}(\lambda)=\frac{1}{\sqrt{2}}\left(\varepsilon_{\lambda_{1}}^{\mu}\varepsilon^{\mu^{\prime}}_{\lambda_{2}}-\varepsilon_{\lambda_{1}}^{\mu^{\prime}}\varepsilon^{\mu}_{\lambda_{2}}\right)\,,$
(67)

and

$(\lambda_{1},\lambda_{2})=(1,0),(1,-1),(0,-1)\,,~{}~{}\text{for}~{}~{}\lambda=1,0,-1\,.$
(68)

To see that $V_{1,2}$ can be viewed as spacelike vectors, in the $\vec{q}$
frame, we define

$\left(\vec{V}_{1}\right)_{i}=V_{1}^{0i}\,,~{}~{}~{}\left(\vec{V}_{2}\right)_{i}=\epsilon_{ijk}\vec{V}_{2}^{jk}\,,$
(69)

with $i,j,k=x,y,z$ and $\epsilon_{ijk}$ the totally anti-symmetric tensor. It
shows that $V_{1}$ and $V_{2}$ are spacelike, which are essentially spin-1
under the $SO(3)$ rotational group.

The results can be understood in terms of the group theory, given as444In the
$SO(3)$ group, $\overline{{\bf 3}}$ and ${\bf 3}$ are equivalent. However,
they behave differently under the $O(3)$ group, where $\overline{{\bf 3}}$ and
${\bf 3}$ are parity even and odd, respectively.

$\left({\bf 1}\oplus{\bf 3}\right)_{1}\otimes\left({\bf 1}\oplus{\bf
3}\right)_{2}=\left({\bf 1}\oplus{\bf 3}\oplus{\bf 6}\right)_{S}+\left({\bf
3}\oplus\overline{{\bf 3}}\right)_{A}\,,$ (70)

where ${\bf 1}$, ${\bf 3}$, and ${\bf 6}$ are the representations of the
$SO(3)$ rotational group, and we have used the fact that a four-vector is
${\bf 1}\oplus{\bf 3}$. In Eq. (70), $``S"$ and $``A"$ in the subscripts
indicate symmetric and antisymmetric between the first and second objects,
respectively. The antisymmetric nature of $\sigma_{\mu\nu}$ forces us to
select the second solution, where ${\bf 3}$ and $\overline{{\bf 3}}$
correspond to $V_{1}$ and $V_{2}$, respectively.

Now, we are able to rewrite the transition matrix element of the tensor
operator as

$\displaystyle
g^{\mu\nu}g^{\mu^{\prime}\nu^{\prime}}\left(\overline{u}_{\ell}\sigma_{\mu\mu^{\prime}}P_{L}v\right)\langle\Lambda_{c}\left|\overline{c}\sigma_{\nu\nu^{\prime}}b\right|\Lambda_{b}\rangle=\sum_{\lambda,n}L^{V_{n}}_{\lambda}B_{\lambda}^{V_{n}}\,,$
$\displaystyle
L^{V_{n}}_{\lambda}=\frac{(-1)^{n}}{2}iV_{n}^{\mu\mu^{\prime}}(\lambda)\overline{u}_{\ell}\sigma_{\mu\mu^{\prime}}P_{L}v\,,\qquad
B_{\lambda}^{V_{n}}=-iV_{n}^{\nu\nu^{\prime}*}(\lambda)\langle\Lambda_{c}\left|\overline{c}\sigma_{\nu\nu^{\prime}}b\right|\Lambda_{b}\rangle\,,$
(71)

which can be understood as a product of $\Lambda_{b}\to\Lambda_{c}V_{n}$ and
$V_{n}\to\ell^{-}\overline{\nu}$, with $V_{n}$ effectively spin-1 particles.
The helicity amplitudes of the lepton sector are given as

$\displaystyle\left(h^{V_{1}}_{+},h^{V_{1}}_{-}\right)=\left(-\frac{1}{\sqrt{2}}h_{-}\,,-\sqrt{2}h_{+}\right)\,,~{}~{}~{}\left(h^{V_{2}}_{+},h^{V_{2}}_{-}\right)=\left(\frac{1}{2}h_{-}\,,h_{+}\right)\,,$
(72)

with

$h_{\pm}^{V_{1,2}}=L^{V_{1,2}}_{\pm\frac{1}{2}-\frac{1}{2}}\left(\lambda_{\ell}=\pm\frac{1}{2},\vec{p}_{\ell}=-\vec{p}_{\nu}=p\hat{z}\right)\,,$
(73)

describing $(V_{1}\to\ell^{-}\overline{\nu})$ and
$(V_{2}\to\ell^{-}\overline{\nu})$, respectively. For the baryon sector, the
helicity amplitudes read as

$H_{\lambda_{c},\lambda}^{V_{n}}=B^{V_{n}}_{\lambda}\left(\lambda_{b}=\lambda_{c}-\lambda,\vec{p}_{c}=-\vec{q}=|\vec{p}_{c}|\hat{z}\right)\,,$
(74)

with $H_{\lambda_{c},t}^{V_{n}}=0$ and $V_{1,2}$ being spacelike.

The six-fold angular distributions now take the form

$\displaystyle\frac{(q^{2}-m_{\ell}^{2})|\vec{p}_{c}|}{q^{2}}$
$\displaystyle\sum_{\lambda_{\ell}\,,\lambda\,,\lambda_{b}}\rho_{\lambda_{b},\lambda_{b}}\left|A^{c}_{\lambda}\sum_{\lambda_{c}\,,\lambda_{W}}\left[C_{T}h_{\lambda_{\ell}}^{V_{1}}\left(H_{\lambda_{c},\lambda_{W}}^{V_{1}}-\frac{1}{\sqrt{2}}H_{\lambda_{c},\lambda_{W}}^{V_{2}}\right)\right.\right.$
$\displaystyle\left.\left.+(-1)^{J_{W}}H_{\lambda_{c}\,,\lambda_{W}}h_{\lambda_{\ell}}\right]d^{\frac{1}{2}}(\theta_{b})^{\lambda_{b}}\,_{\lambda_{c}-\lambda_{W}}d^{\frac{1}{2}}(\theta_{c})^{\lambda_{c}}\,_{\lambda}d^{J_{W}}(\theta_{\ell})^{\lambda_{W}}\,_{\lambda_{\ell}-\frac{1}{2}}e^{i(\lambda_{c}\phi_{c}+\lambda_{\ell}\phi_{\ell})}\right|^{2}\,,$

which can not be reduced to Eq. (II) by redefining the amplitudes. Thus, Table
1 would no longer be suitable after the tensor operator is considered.

## References

  * (1) M. Jung and D. M. Straub, JHEP 01, 009 (2019); Q. Y. Hu, X. Q. Li, X. L. Mu, Y. D. Yang and D. H. Zheng, JHEP 06, 075 (2021); Z. R. Huang, E. Kou, C. D. Lü and R. Y. Tang, Phys. Rev. D 105, 013010 (2022); K. Ezzat, G. Faisel and S. Khalil, arXiv:2204.10922; A. Datta, H. Liu and D. Marfatia, arXiv:2204.01818; R. Y. Tang, Z. R. Huang, C. D. Lü and R. Zhu, arXiv:2204.04357.
  * (2) R. Aaij et al. [LHCb], Phys. Rev. Lett. 115, 111803 (2015) [erratum: Phys. Rev. Lett. 115, 159901; S. Hirose et al. [Belle], Phys. Rev. Lett. 118, 211801 (2017); R. Aaij et al. [LHCb], Phys. Rev. Lett. 120, 171802 (2018).
  * (3) P. A. Zyla et al. [Particle Data Group], PTEP 2020, 083C01 (2020).
  * (4) Y. S. Amhis et al. [HFLAV], Eur. Phys. J. C 81, 226 (2021).
  * (5) S. Jaiswal, S. Nandi and S. K. Patra, JHEP 12, 060 (2017).
  * (6) S. Bifani, S. Descotes-Genon, A. Romero Vidal and M. H. Schune, J. Phys. G 46, 023001 (2019).
  * (7) R. Aaij et al. [LHCb], Phys. Rev. Lett. 128, 191803 (2022).
  * (8) F. U. Bernlochner, Z. Ligeti, D. J. Robinson and W. L. Sutcliffe, Phys. Rev. Lett. 121, 202001 (2018); F. U. Bernlochner, Z. Ligeti, D. J. Robinson and W. L. Sutcliffe, Phys. Rev. D 99, 055008 (2019).
  * (9) S. Shivashankara, W. Wu and A. Datta, Phys. Rev. D 91, 115003 (2015); R. Dutta, Phys. Rev. D 93, 054003 (2016); X. Q. Li, Y. D. Yang and X. Zhang, JHEP 02, 068 (2017); A. Datta, S. Kamali, S. Meinel and A. Rashed, JHEP 08, 131 (2017); E. Di Salvo, F. Fontanelli and Z. J. Ajaltouni, Int. J. Mod. Phys. A 33, 1850169 (2018); X. L. Mu, Y. Li, Z. T. Zou and B. Zhu, Phys. Rev. D 100, 113004 (2019); P. Böer, A. Kokulu, J. N. Toelstede and D. van Dyk, JHEP 12, 082 (2019); X. L. Mu, Y. Li, Z. T. Zou and B. Zhu, Phys. Rev. D 100, 113004 (2019); N. Penalva, E. Hernández and J. Nieves, Phys. Rev. D 100, 113007 (2019); P. Colangelo, F. De Fazio and F. Loparco, JHEP 11, 032 (2020); N. Penalva, E. Hernández and J. Nieves, JHEP 06, 118 (2021); N. Penalva, E. Hernández and J. Nieves, Phys. Rev. D 101, 113004 (2020).
  * (10) A. Kadeer, J. G. Korner and U. Moosbrugger, Eur. Phys. J. C 59, 27 (2009); T. Gutsche, M. A. Ivanov, J. G. Körner, V. E. Lyubovitskij, P. Santorelli and N. Habyl, Phys. Rev. D 91, 074001 (2015).
  * (11) Q. Y. Hu, X. Q. Li, Y. D. Yang and D. H. Zheng, JHEP 02, 183 (2021); N. Penalva, E. Hernández and J. Nieves, JHEP 04, 026 (2022).
  * (12) R. Aaij et al. [LHCb], JHEP 06, 110 (2020).
  * (13) M. Ferrillo, A. Mathad, P. Owen and N. Serra, JHEP 12, 148 (2019).
  * (14) S. Iguro and R. Watanabe, JHEP 08, 006 (2020); S. Iguro, M. Takeuchi and R. Watanabe, Eur. Phys. J. C 81, 406 (2021).
  * (15) S. Groote, J. G. Körner and B. Melić, Eur. Phys. J. C 79, 948 (2019).
  * (16) J. Y. Cen, C. Q. Geng, C. W. Liu and T. H. Tsai, Eur. Phys. J. C 79, 946 (2019).
  * (17) C. Q. Geng, C. W. Liu and T. H. Tsai, Phys. Rev. D 103, 054018 (2021).
  * (18) C. Q. Geng and C. W. Liu, JHEP 11, 104 (2021).
  * (19) C. W. Liu and C. Q. Geng, arXiv:2205.08158.


# Constrained simulations of the local Universe with Modified Gravity

Krishna Naidoo<EMAIL_ADDRESS>Center for Theoretical Physics, Polish
Academy of Sciences, al. Lotników 32/46 Warsaw, Poland  Department of Physics
& Astronomy, University College London, Gower Street, London, WC1E 6BT, UK
Wojciech A. Hellwing<EMAIL_ADDRESS>Center for Theoretical Physics,
Polish Academy of Sciences, al. Lotników 32/46 Warsaw, Poland  Maciej Bilicki
Center for Theoretical Physics, Polish Academy of Sciences, al. Lotników 32/46
Warsaw, Poland  Noam Libeskind  Leibniz-Institut für Astrophysik Potsdam, An
der Sternwarte 16, D-14482 Potsdam, Germany  University of Lyon, UCB Lyon 1,
CNRS/IN2P3, IUF, IP2I Lyon, France  Simon Pfeifer  Leibniz-Institut für
Astrophysik Potsdam, An der Sternwarte 16, D-14482 Potsdam, Germany  Yehuda
Hoffman  Racah Institute of Physics, Hebrew University, Jerusalem 91904,
Israel

###### Abstract

We present a methodology for constructing modified gravity (MG) constrained
simulations of the local Universe using positions and peculiar velocities from
the CosmicFlows data set. Our analysis focuses on the following MG models: the
normal branch of the Dvali-Gabadadze-Porrati (nDGP) model and Hu-Sawicki
$f(R)$ model. We develop a model independent methodology for constructing
constrained simulations with any given power spectra and numerically
calculated linear growth functions. Initial conditions (ICs) for a set of
constrained simulations are constructed for the standard cosmological model
$\Lambda$CDM and the MG models. Differences between the model’s reconstructed
Wiener filtered density and the resultant simulation density are presented
showing the importance for the generation of MG constrained ICs to study the
subtle effects of MG in the local Universe. These are the first MG constrained
simulations ever produced. The current work paves the way to improved
approximate methods for models with scale-dependent growth functions, such as
$f(R)$, and for high-resolution hydrodynamical MG zoom-in simulations of the
local Universe.

##  I Introduction

The standard model of cosmology, $\Lambda$ \+ cold dark matter ($\Lambda$CDM)
stands on firm pillars formed by a multitude of observational tests. From the
early epochs of primordial nucleosynthesis [1], and the statistical properties
of the cosmic microwave background radiation and its angular fluctuations [2],
to the late time formation and evolution of large scale structure as
gloriously manifested by the amazing spatial patterns observed from vast and
deep modern-age galaxy surveys (such as the Baryon Oscillation Spectroscopic
Survey [3]). The observations concerning all of these epochs and phenomena can
be explained remarkably well by a simple 6-parameter $\Lambda$CDM model. These
undeniable empirical successes that form the foundations of $\Lambda$CDM come
however mostly either from the large scales or early epochs when the physics
of the involved phenomena is generally in the linear or mildly-nonlinear
regimes.

In the past few decades the precision, quality and volume of the data
concerning the so-called local Universe has grown and improved many-folds.
Thus opening the window on non-linear scales, where the environment defined by
the large-scale matter/galaxy distribution is entangled with local non-linear
processes driving the evolution of galaxies, and their motions and clustering.
Dedicated observational campaigns have brought us a detailed and interesting
picture of our close cosmic neighbourhood – within a region of 100-200
$h^{-1}{\rm Mpc}$ from the Local Group. With new observations of hundreds and
thousands of local dwarf galaxies, more precise measurements of the Virgo,
Coma and other local galaxy clusters to the stunning Cosmic Flows data
containing distances and velocities of tens of thousands of nearby galaxies.
All these new local Universe data have a great potential for providing new
stringent tests for cosmology and new insights into non-linear phenomena and
physics of galaxy formation and evolution. They also provide new intricate
details concerning the local Universe that require robust modelling and
understanding of structure formation physics in the non-linear regime.

Classically, the high-resolution state-of-the-art cosmological simulations
provide a powerful tool for modeling and understanding the physics of the non-
linear structure formation regime [4]. However, since classical simulations
use random phases of the initial conditions, the structures reflect just one
particular random realisation from a vast ocean of potential configurations.
Aiming at simulating the very local universe requires running very high
resolution simulations, to capture the intricate physics in the deeply non-
linear regime. At the same time one needs to run large box simulations, in
order to sample Local Group-like environments out of cosmic variance. Large
box high resolution galaxy formation simulations are prohibitively expensive.

Here, the constrained simulations [5] of the local Universe which reproduce
structures and phase information from observations, allows us to circumvent
the aforementioned limitations of random-phase simulations. These new class of
carefully designed and engineered simulations nowadays play a critical role in
distinguishing the subtle features from model extensions to $\Lambda$CDM. In
constrained simulations the variance of the structures and resultant
statistics is greatly reduced – a property which is especially important in
the context of the local Universe where observations of galaxies are made to
higher accuracy (than those made on distant galaxies typically measured from
galaxy redshift surveys) and there is greater sensitivity to small scale
physics (such as the properties of dwarf galaxies) that is not available to
larger surveys. Simulations mimicking structures from real galaxy observations
have been produced from a variety of methodologies; constrained simulations
[5] using peculiar velocities constraints from the local Universe [6, 7, 8]
(the focus of this paper), Bayesian reconstruction of the density field [9,
10] and methods using both density and peculiar velocity constraints [11].

In this paper we focus on constructing constrained simulations of the local
Universe using peculiar velocities [12]. Unlike estimations of the density
field from galaxy positions – biased tracers of the density field – peculiar
velocities in the linear regime are related directly to the density field.
This allows us to estimate the density field with limited dependence on the
galaxy tracer bias, while the bias from non-linear and stochastic galaxy
motions are limited by considering only large scale linear regime
reconstructions of the density field. Furthermore, peculiar velocities come
with the following limitations: (1) measurements are limited to only the
radial component of the velocities, (2) they can only be measured at low
redshift (depending on accurate distance measurements and assumptions on the
background expansion and Hubble constant) (4) non-linear velocities need to be
removed and (5) Malmquist bias need to be considered. See [13] for a more
detailed review. As a result constrained simulations from peculiar velocities
have been limited to the local environment, in particular using data from
CosmicFlows [see 14]. Several studies, such as the Constrained Local UniversE
Simulations (CLUES) [6, 7] and the High-resolution Environmental Simulations
of The Immediate Area (HESTIA) [8], have used peculiar velocities from
CosmicFlows to produce constrained simulations of the local Universe. However,
these simulations have been limited to $\Lambda$CDM and the model dependence
of the reconstruction and the reproduced structures has yet to be explored.

Modified gravity (MG) are a class of models that extend Einstein’s theory of
general relativity. They are motivated by a desire to provide a better
theoretical explanation of the cosmological constant, a constituent of nature
that has proven to be a significant challenge for theoretical fundamental
physics. They achieve this by introducing an effective fifth force, that
replicates the effects of the cosmological constant but that is screened on
relatively small cosmological scales. The range of this screening mechanism
and the high fidelity observations of the local Universe, make it a perfect
environment for extracting the subtle effects of MG. In this paper we focus on
the following MG models: Hu-Sawicki $f(R)$ with chameleon screening [15] and
the normal branch of the Dvali-Gabadadze-Porrati (nDGP) with Vainshtein
screening [16]. Constrained simulations have never been produced for models of
MG and those produced for the local Universe would enable detailed comparison
studies looking to smaller scales than is typically accessible to large
cosmological galaxy surveys.

In this paper we extend the methodology used by CLUES for creating constrained
local Universe simulations to two models of MG. We particularly focus on
making these methods numerically driven, enabling the extension to other
models in future studies. The methods used by CLUES follows the procedures in
[17, 18, 19] which are implemented in ICeCoRe (Initial Conditions &
Constrained Realisations) which is used to generate initial conditions (ICs).
To construct ICs we first use the Wiener filtered [WF; 20] to reconstruct the
density field which is used to calculate the linear displacement field. The
displacement field is used to push the constraint points to their locations in
Lagrangian space at early times. This assumes linear perturbation theory and
is referred to as the reverse Zeldovich approximation [RZA; 17]. Constrained
realisations are then generated [5, 19] from the corrected location of
galaxies after applying the RZA and the density scaled to an initial redshift
to which ICs are generated. In this work we modify these methods to allow for
MG ICs. We test and run constrained simulations using the COmoving Lagrangian
Acceleration method [COLA; 21], specifically MG-PICOLA [22] which implements
the COLA method for a variety of MG models.

The paper is organised as follows, in Sec. II we discuss the theory behinds
the models and constrained realisations, in Sec. III we discuss the
methodology for the current $\Lambda$CDM implementation and modification for
MG, in Sec. IV we discuss the results on the constrained initial conditions
and constrained COLA simulations and in Sec. V we discuss the results, future
work and challenges.

##  II Theory

###  II.1 Cosmological models

Here, we provide a brief description of the cosmological models we
investigate. Our fiducial case that we take as a baseline with respect to
which we measure all deviations and signals is the standard $\Lambda$CDM. We
explore the physics of non-standard structure formation in two families of MG
models, $f(R)$ and nDGP. For all the simulations we assume the same background
cosmology based on [2] with parameter values: $\Omega_{\rm m}=0.3111$,
$\Omega_{\rm b}=0.049$, $\Omega_{\Lambda}=0.6889$, $H_{0}=67.66$, $A_{\rm
s}=2.105\times 10^{-9}$ and $n_{\rm s}=0.9665$ in a spatially flat universe
(_i.e._ $\Omega_{\rm k}=0$) and with neutrinos assumed to be massless (_i.e._
$\sum m_{\rm\nu}=0\,{\rm eV}$). That is all models share the same expansion
history, and the various cosmologies differ in the structure formation pace
(_i.e._ the linear growth rate) and the locally operating highly non-linear
fifth-forces and screening mechanism (for the case of MG).

####  II.1.1 Modified gravity models: nDGP and f(R)

Observational constraints from the Solar System and massive bodies
significantly reduce the possible extensions and deviations from Einstein’s
theory of General Relativity. Any modifications need to simultaneously match
Solar System constraints for gravity but at the same time need to allow for
possible departures from $\Lambda$CDM on large cosmological scales. Typically
this is achieved in MG by the inclusion of a screening mechanism.

In nDGP, gravity is propagated through extra dimensions, unlike other forces,
while in $f(R)$, non-linear functions are added to the Ricci scalar. In both
cases, the new effective action integral will allow for extra degrees of
freedom, that can be modelled by additional scalar fields and their
interactions with matter. As a result, the new dynamics of these class of
models allows for a non-vanishing fifth scalar-like force operating on
cosmological and intergalactic scales. The physics of models from both
families naturally contains a non-linear mechanism to suppress propagation of
such fifth-forces. These are called screening mechanism, and in general the
effective range of their operation is limited to small non-linear scales. In
the nDGP model, screening is achieved by the means of the Vainshtein effect
[23]; while in the $f(R)$ family, the corresponding mechanism is the Chameleon
screening effect [24]. The Vainshtein screening is dependent on the mass and
distance from a screened object with no explicit dependence on either local or
global environment, while the chameleon screening is dependent on the local
curvature, and thus effectively on the density distribution of the local
matter fields.

As a test-bed for the MG theories we select one particular model for nDGP and
one model for $f(R)$. In both cases we assume a background cosmology following
the parameters given for the fiducial $\Lambda$CDM model discussed above with
the additional parameter $r_{\rm c}H_{0}=1$ for nDGP and $|f_{R0}|=10^{-5}$
for $f(R)$. We refer to these models as N1 and F5, respectively. At high
redshift these models exhibit almost identical clustering properties as
$\Lambda$CDM but depart at low redshift. This is showcased by the differences
in their linear power spectra at redshift $z=0$ (shown in Fig. 1 and obtained
from a modified version of CAMB) where nDGP shows a scale independent shift,
while in $f(R)$ the shift is scale dependent – with large scales consistent
with $\Lambda$CDM but departing at smaller scales.

Figure 1: The linear power spectra at redshift $z=0$ is shown for the fiducial
$\Lambda$CDM model, the nDGP model N1, and the $f(R)$ model F5 in the top
panel. Differences with respect to $\Lambda$CDM are shown in the bottom panel.

###  II.2 Constrained simulations

The motivations for constrained simulations are many. One particular goal is
to create simulations that mimic properties of data and the view of our local
Universe that it offers. Another important aspect of constrained simulations
is that they allow for studying, and thus limiting to some extent, the impact
of cosmic variance on nearby non-linear structures. In this context, it is
crucial to remember that the observational galaxy data is noisy and sparse.
Thus the constrained simulations both replicate properties of the data but
also add information where data is missing.

Our approach involves three major steps: (1) WF which reconstructs the linear
velocity and density field from peculiar velocity measurements; (2) the RZA
which reduces linear-order shifts in the position of constraints relative to
an initial redshift, and (3) the generation of constrained realisations which
adds random fluctuations in places which are poorly constrained. It is
important to note that constrained realisations assume the input data and
realisations are Gaussian. Of course low redshift observations of large scale
structure are certainly non-Gaussian, so some care is needed to try to limit
measurements to the linear regime. The latter is usually attained by
considering all the relevant density and velocity fields smoothed on scales
where the effects of the non-linearities are already reduced. The procedure
for obtaining constrained field realisations using the observed velocity data
is described in detail in [12]. Below, for completeness, we discuss the main
steps of this method.

Figure 2: The density and velocity correlation functions are shown for the
fiducial $\Lambda$CDM model and the MG models of nDGP N1 and $f(R)$ F5. In the
top left panel $\xi$ is the density-density correlation function, in the top
right panel $\zeta$ is the velocity-density cross correlation function and in
the bottom panels are the velocity-velocity correlation functions ($\psi_{\rm
R}$ the radial component on the left and $\psi_{\rm T}$ the tangential
component on the right). Subpanels show the differences between these
correlation functions with respect to the base fiducial $\Lambda$CDM model.

####  II.2.1 Wiener filtering

Following [12] the WF density $\delta^{\rm WF}$ is estimated from a given set
of peculiar velocity $\boldsymbol{u}$ constraints from the relation

$\delta^{\rm
WF}(\boldsymbol{r})=\sum_{i=1}^{M}\sum_{j=1}^{M}\langle\delta(\boldsymbol{r})\,\boldsymbol{u}_{i}\rangle\,\langle\boldsymbol{u}_{i}\,\boldsymbol{u}_{j}\rangle^{-1}\,\boldsymbol{u}_{j}\,,$
(1)

where $\langle\delta(\boldsymbol{r})\,\boldsymbol{u}_{i}\rangle$ is given by
the linear theory velocity-density correlation,
$\langle\boldsymbol{u}_{i}\,\boldsymbol{u}_{j}\rangle$ is the peculiar
velocity covariance matrix and the indices $i$ and $j$ sum over $M$ peculiar
velocity constraints. The density-velocity correlation vector is defined as

$\langle\delta(\boldsymbol{r})\,\boldsymbol{u}_{i}\rangle=-\dot{a}f\zeta(r^{\prime})\left(\boldsymbol{e}_{i}\cdot\boldsymbol{\bar{r}}^{\prime}\right),$
(2)

where $a$ is the scale factor related to the redshift $z$ by $a=1/(1+z)$,
$\dot{a}=a\,H(z)$ and $H(z)$ is the Hubble parameter (note at $z=0$ this
reduces to $\dot{a}=H_{0}$). The logarithmic growth rate $f={\rm d}\ln D/{\rm
d}\ln a$ where $D$ is the linear growth function normalised at $z=0$ and
$\zeta(r)$ is the density-velocity correlation function. Lastly
$r^{\prime}=|\boldsymbol{r}^{\prime}|$,
$\boldsymbol{\bar{r}}^{\prime}=\boldsymbol{r}^{\prime}/r^{\prime}$ and
$\boldsymbol{e}_{i}$ is a unit vector describing the direction of the peculiar
velocity constraint $\boldsymbol{u}_{i}$, in this study this is equivalent to
the line-of-sight vector.

The velocity covariance matrix is given by

$\langle\boldsymbol{u}_{i}\,\boldsymbol{u}_{j}\rangle=\boldsymbol{e}_{i}\cdot\langle\boldsymbol{u}_{i}(\boldsymbol{r})\,\boldsymbol{u}_{j}(\boldsymbol{r}+\boldsymbol{r}^{\prime})\rangle\cdot\boldsymbol{e}_{j}+\sigma_{i}^{2}\delta^{\rm
K}_{ij},$ (3)

where $\sigma_{i}$ are measurement errors on the peculiar velocities,
$\delta^{\rm K}_{ij}$ is the Kronecker delta function and the the components
of the linear theory peculiar velocity covariance matrix (_i.e._ velocity
correlation tensor) are given by

$\begin{split}\Psi^{v}_{\alpha\beta}\equiv\langle\boldsymbol{u}(\boldsymbol{r})\,&\boldsymbol{u}(\boldsymbol{r}+\boldsymbol{r}^{\prime})\rangle_{\alpha\beta}\,=(\dot{a}f)^{2}\\\
&\times\bigg{\\{}\psi_{\perp}(r^{\prime})\delta^{K}_{\alpha\beta}+\big{[}\psi_{\parallel}(r^{\prime})-\psi_{\perp}(r^{\prime})\big{]}\bar{r}^{\prime}_{\alpha}\bar{r}^{\prime}_{\beta}\bigg{\\}}.\end{split}$
(4)

Where $\psi_{\perp}(r)$ and $\psi_{\parallel}(r)$ are the tangential and
radial velocity-velocity correlation functions. The above decomposition of the
full velocity correlation into its radial and tangential holds for a
statistically homogeneous and isotropic velocity field [25]. We can define the
auto- and cross-correlation functions for the density and peculiar velocities
as follows

$\displaystyle\xi(r)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}k^{2}\,P(k)\,j_{0}(kr)\,{\rm
d}k\,,$ (5)
$\displaystyle\zeta(r)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}k\,P(k)\,j_{1}(kr)\,{\rm
d}k\,,$ (6)
$\displaystyle\psi_{\parallel}(r)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\left(P(k)\,j_{0}(kr)-2\,\frac{j_{1}(kr)}{kr}\right)\,{\rm
d}k\,,$ (7)
$\displaystyle\psi_{\perp}(r)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}P(k)\,\frac{j_{1}(kr)}{kr}\,{\rm
d}k\,,$ (8)

where $\xi(r)$ is the density-density correlation function. The integrals are
functions of the linear matter power spectra $P(k)$ and spherical Bessel
functions $j_{n}(x)$ (_i.e._ a spherical Bessel function at $x$ of order $n$).
In the above relation, we will be using consistently $P(k)$ and $f$ of either
$\Lambda$CDM or a given MG model respectively. The numerically computed
correlation functions are shown in Fig. 2 for two MG models and $\Lambda$CDM.
Negative correlations are seen for large scales in $\xi$ and $\psi_{\rm R}$
and the wiggles seen in the subpanel comparisons of $\xi$ originate from
numerical instabilities at larger radii where the correlation values are
small.

####  II.2.2 Reverse Zeldovich approximation

If constrained simulations were generated directly from the peculiar velocity
constraints at their native (_i.e._ final observed) positions, the simulations
generated would have structures that are systematically shifted from the input
field. This is because the location of structures move with respect to an
arbitrary point of origin as they grow and evolve. This is a direct
consequence of large and small-scale coherent (_i.e._ bulk) flow motions in
the Universe. In order to account for this and to ensure that the final
simulated structures appear at the comoving positions corresponding with the
observational data locations we need to push our constraints to where they
would appear at early times _i.e._ the initial redshift. In [17] the dominant
effect was shown to be well captured already by the terms from linear
perturbation theory. Thus it is sufficient to perform a simple linear order
shift to the positions of the peculiar velocity constraints, known as the
Reverse Zeldovich approximation (RZA) [26].

The Zeldovich approximation [27], based on linear-order perturbation theory,
can be used to describe the relation between the initial (_i.e._ Lagrangian)
position of a fluid element in space, with its final (_i.e._ Eulerian)
position as a function of time

$\boldsymbol{x}(\boldsymbol{q},z)=\boldsymbol{q}+\boldsymbol{\psi}(\boldsymbol{q},z)\,,$
(9)

where $\boldsymbol{x}$ is the final position in space, $\boldsymbol{q}$ the
initial condition position and $\boldsymbol{\psi}$ the displacement field, all
given in comoving coordinates. The peculiar velocities $\boldsymbol{u}$ are
related to the displacement field by

$\boldsymbol{u}(\boldsymbol{q},z)=\dot{a}\,f\,\boldsymbol{\psi}(\boldsymbol{q},z)\,.$
(10)

In the RZA the following approximation is employed,

$\boldsymbol{q}_{\rm
RZA}(\boldsymbol{x})\sim\boldsymbol{x}(z)-{\boldsymbol{v}(\boldsymbol{x},z)\over\dot{a}f}\,,$
(11)

to give the initial condition positions of the peculiar velocity constraints.
Note, in our work $\boldsymbol{\psi}$ is the linear displacement computed from
the reconstructed WF velocity field.

####  II.2.3 Constrained Realisations

The final step is to now construct specific initial condition realisations
from the constraints (_i.e._ the data) after performing the RZA. This creates
realisations which retain information at locations with high constraints, but
with added noise and random information at scales and locations which are
poorly constrained. The density contrast for a constrained realisation (CR) is
obtained from the following relation

$\begin{split}\delta^{\rm CR}(\boldsymbol{r})&=\delta^{\rm
RR}(\boldsymbol{r})\\\
&+\sum_{i=1}^{M}\sum_{j=1}^{M}\langle\delta(\boldsymbol{r})\,\boldsymbol{u}_{i}\rangle\,\langle\boldsymbol{u}_{i}\,\boldsymbol{u}_{j}\rangle^{-1}\,(\boldsymbol{u}_{j}-\boldsymbol{u}^{\rm
RR}_{j})\,,\end{split}$ (12)

following [5] where $\boldsymbol{u}$ are the peculiar velocities at their
initial condition positions (i.e. after applying the RZA), $\delta^{\rm RR}$
the density from a random realisation and $\boldsymbol{u}^{\rm RR}$ the
peculiar velocity from the random realisation at the position and direction of
the peculiar velocity constraints. The density and displacement field are
related in Fourier space by

$\boldsymbol{\psi}(\boldsymbol{k})=\frac{i\boldsymbol{k}}{k^{2}}\delta(\boldsymbol{k}),$
(13)

since $\nabla^{2}\Phi=\delta$ and
$\boldsymbol{\psi}=-\boldsymbol{\nabla}\Phi$, and where $k=|\boldsymbol{k}|$
and $i$ is the imaginary unit.

####  II.2.4 ICeCoRe $f(R)$ implementation

The production of constrained realisations from ICeCoRe broadly follows the
steps outlined above. However, ICeCoRe was designed for constructing
constrained realisations in $\Lambda$CDM and as a result there are a few
design choices which are $\Lambda$CDM-specific. Although these limitations do
not limit the accuracy of constrained nDGP simulations they do affect the
accuracy of constrained $f(R)$ simulations. The most important of them is the
choice to take constraints of the displacement field rather than the velocity
field directly. The two are linked by Eq. (10) where a simple division by the
factor $\dot{a}f$ can convert $\boldsymbol{u}$ to $\boldsymbol{\psi}$, but
this is only true when $f$ is scale independent. To correctly convert one to
the other in $f(R)$ would require a multiplication in Fourier space, but this
is impractical given the constraints $\boldsymbol{u}$ do not lie on a uniform
grid (where we can take advantage of Fast Fourier Transforms – FFT and
numerical methods). The choice to constrain the displacement field also leads
to the simplification of Eq. (2) and Eq. (4) where the $\dot{a}f$ terms drop
from the equation. To generate self-consistent $f(R)$ constrained simulations
will require developing ICeCoRe with a different philosophy. This will require
moving the linear growth function $f$ from Eqs. 2 and 4 to inside the
integrals of Eqs. (6–8). Furthermore, we will need to use the peculiar
velocity constraints directly to generate the constrained realisations rather
than converting to the displacement field. Implementing this change will be
the subject of future work. For the moment we instead consider an effective
scalar value for $f$ at some chosen scale.

##  III Method

In this section we discuss the steps for constructing constrained simulations
for $\Lambda$CDM and MG. We begin by discussing the peculiar velocity
constraints used, then discuss the implementation of ICeCoRe for generating
initial conditions for $\Lambda$CDM followed by a discussion of the
modifications to the ICeCoRe inputs and intermediate procedures for generating
initial conditions for MG. Lastly, we discuss the subsequent COLA $N$-body
simulations and the random realisations used for comparisons.

###  III.1 Peculiar velocity constraints from CosmicFlows-2

We use CosmicFlows-2 (CF2) [28], a catalog of galaxy redshifts and distances,
and hence of galaxy distances and peculiar velocities, for imposing the
constraints. The data consist of measurements of the sky positions, distance
moduli and redshifts for over 10,000 local Universe galaxies and groups of
galaxies, which are mostly contained within a distance of $\lesssim
150\,h^{-1}{\rm Mpc}$. Galaxies which are closely clustered are grouped and
the peculiar velocity constraints replaced by their mean – suppressing virial
motions in high density environments [29] and maintaining the suitability of
linear theory assumptions. A straight forward transformation from the noisy
distance muduli and redshifts to distances and velocities induces Malmquist-
like biases. The data used here is corrected by the bias minimization
algorithm [30].

Although newer CosmicFlows data releases are available [31, 14], CF2 [28] has
been widely used in the community for the production of constrained
simulations by CLUES [7] and systematics are well understood [30], while the
effects in later releases still require further study. As we are simply
interested in comparing and studying the implementation of MG for constrained
simulations the most up-to-date data are not required for the time being.

###  III.2 Constrained realisations for $\Lambda$CDM

Constrained simulations are generated by ICeCoRe in four steps. In each step
ICeCoRe takes as input the linear power spectrum $P(k)$ which is computed for
$\Lambda$CDM using the cosmological Boltzmann code solver CAMB
111https://camb.readthedocs.io/.

  1. 1.

Wiener Filtering: Peculiar velocity constraints from CF2 are used to construct
the WF density field and velocity field computed by differentiating the WF
density field in Fourier space.

  2. 2.

Reverse Zeldovich Approximation: Using interpolated values of the displacement
field, from the WF reconstruction, we subtract the linear-order displacement
to place the constraints to their positions in the initial conditions in
Lagrangian space.

  3. 3.

Constrained realisations: A random seed is given to generate a CR of the
density field using the peculiar velocities at the RZA positions following the
density power spectrum at $z=0$.

  4. 4.

Initial conditions: The density field at the initial conditions are calculated
by multiplying the CR density by the linear growth function at the initial
redshift $z_{\rm IC}=49$. IC particles are then generated on a grid using the
Zeldovich approximation.

Figure 3: The scale-independent growth functions $D(z)$ and $f(z)$ are shown
for the fiducial base $\Lambda$CDM model (in black), the nDGP model N1 (in
blue) and the scale-dependent growth functions shown for the $f(R)$ model F5
(where colours correspond to scales $k$ in the colourbar) and the red line
shows the effective growth functions at a scale of $k=2\pi/12\,{h{\rm
Mpc}^{-1}}$.

###  III.3 Constrained realisations for modified gravity

Constrained simulations for MG are implemented with ICeCoRe in a similar
method to $\Lambda$CDM. We take advantage of the modular set up to alter the
inputs and intermediate steps for MG and compute scale independent (for nDGP)
and dependent (for $f(R)$) growth functions numerically from MG power spectra.

The first modification is the input of a MG power spectra. These are computed
using MGCAMB 222Made available to us by Hans Winter a modified version of CAMB
used to compute power spectra for a range of MG models. The second step
requires the pre-calculation of growth functions. Rather than hard-code this
procedure into our pipelines, we compute growth functions completely
numerically from linear power spectra evaluated at many redshifts. The linear
scale-dependent growth function is calculated by taking

$D(k,z)=\sqrt{\frac{P(k,0)}{P(k,z)}},$ (14)

and the logarithmic growth function computed by numerical differentiation,

$f(k,z)=\frac{{\rm d}\log D(k,z)}{{\rm d}\log a}.$ (15)

Note, for the above we first interpolate $D(k,z)$ onto a regular $\log a$ grid
using a cubic spline, since numerical differentiation performs best on a
regular grid. For the scale-independent case we take the mean of $k$ above
superhorizon scales, _i.e._ $0.01<k<10\,h{\rm Mpc}^{-1}$ (since we are using
the linear power spectra the inclusion of non-linear scales has no effect on
the linear growth functions except to provide more stable and accurate
results). In Fig. 3 we show the numerically calculated growth functions for
nDGP (scale-independent) and $f(R)$ (scale-dependent).

Figure 4: The Wiener filtered (WF) reconstruction of the density field from
peculiar velocity data from CosmicFlows-2 is shown for $\Lambda$CDM on the
bottom left, the N1 nDGP model on the bottom middle and the F5 $f(R)$ model on
the bottom right panel. In the top panels we show deviations of the WF
reconstruction of N1 (middle) and F5 (right) from $\Lambda$CDM. The WF
recovers the linear density field extrapolated to the present epoch, hence the
unphysical values below $-1$. The maps are shown with respect to the
supergalactic $x$ (SGX) and $y$ (SGY) coordinate axis for a slice of width
$20\,{h^{-1}{\rm Mpc}}$ in the supergalactic $z$ axis (_i.e._ $-10\leq{\rm
SGZ}\leq 10\,{h^{-1}{\rm Mpc}}$). The WF reconstruction show fairly consistent
structures with the strongest deviations seen from F5.

The steps to produce constrained realisations for MG follow the procedure
outlined for $\Lambda$CDM with some important differences: such as the use of
MG power spectra and numerically computed growth functions. Note, for $f(R)$
we input growth functions at a single scale, $k=2\pi/12\,h{\rm Mpc}^{-1}$,
since scale dependent growth functions are currently not implemented in
ICeCoRe. This choice of scale is motivated by the effective smoothing of the
Wiener filter reconstruction (of the density field) from peculiar velocities
which is roughly on the scale of $12\,h^{-1}{\rm Mpc}$. The other additional
change we make is to rescale the CR densities to the IC redshift using scale-
dependent growth functions via FFT, which is implemented outside of ICeCoRe.

Figure 5: The density of constrained simulations for $\Lambda$CDM and MG are
shown in a region of 100 $h^{-1}{\rm Mpc}$ around the observer. On the bottom
panels are constrained simulations all generated with $\Lambda$CDM ICs while
on the top are generated with MG ICs. The difference between the structures
from different models is hard to distinguish by eye; for visual comparisons we
indicate the differences between positions of matched particles with respect
to $\Lambda$CDM with arrows/quivers.

###  III.4 Constrained and random simulations

To test the effect of applying constraints for different models we also
construct initial conditions with the same random seeds with the corresponding
$\Lambda$CDM or MG power spectra. Once these are produced we then use the
initial condition particles to run COLA simulations using MG-PICOLA which
implements the COLA method for the MG models of nDGP and $f(R)$ [22]. To be
able to construct halo catalogues from COLA requires that the particle-mesh
grid be a factor of $>3$ times the number of particles along a single axis.
For most of the simulations our initial conditions are generated with
$512^{3}$ particles over a $500\,h^{-1}{\rm Mpc}$ box, so for added validity
and due to the relative inexpensive computation cost of running COLA we use a
particle-mesh grid size of $2048^{3}$. For each model we produce 5 constrained
initial conditions and 5 random initial conditions with the same random seeds.
Constrained simulations are then computed with MG-PICOLA for each model. For
MG models these are run from both $\Lambda$CDM and MG initial conditions. In
total this means the construction of 30 initial conditions from which 50
simulations were produced.

##  IV Results

In this section we compare results obtained from constrained initial
conditions from $\Lambda$CDM and MG models nDGP and $f(R)$; including the WF
reconstruction and COLA simulation outputs and statistics.

###  IV.1 Comparing Wiener filtered reconstruction

The WF reconstruction of the density field computed by ICeCoRe is shown in
Fig. 4. We see the fields are strikingly similar, with the differences
originating from differences in the model’s power spectra and correlation
functions (shown in Fig. 1 and 2, respectively). For nDGP the variations
appear equal in amplitude throughout the volume, perhaps unsurprising given
the differences between the power spectra of $\Lambda$CDM can be described by
an effective amplitude shift on all scales. For $f(R)$ the differences to
$\Lambda$CDM are strongest closer to the observer, showing a strong relation
to the amplitude of the density field.

###  IV.2 Comparing COLA simulation statistics

Figure 6: The density (top) and radial velocity (bottom) fields are compared
between MG and $\Lambda$CDM constrained simulations with MG ICs (left) and
$\Lambda$CDM ICs (right). The density and velocity fields are computed on a
$256^{3}$ grid using the triangular-shaped-cloud mass-assignment scheme and
compared within a radius of $30\,h^{-1}{\rm Mpc}$ from the origin _i.e._ where
most of the CF2 constraints are located. On the y-axis we show the difference
in densities between $\Lambda$CDM and MG on the top row and radial velocities
on the bottom row.

To test the validity of the constrained simulation methodology for MG we began
by comparing the power spectra and halo mass functions computed from the
simulations. These were found to be completely consistent between constrained
and random simulations. Since these statistics are insensitive to phases, this
result is not surprising and we move onto comparisons of the density and
velocity fields.

In Fig. 5 we compare the constrained simulation density fields for the first
realisation. Particles within $\pm 10\,h^{-1}{\rm Mpc}$ from ${\rm SGZ}=0$ at
redshift $z=0$ are used to construct a 2-dimensional projection of the density
on a slice of width $100\,h^{-1}{\rm Mpc}$ around the origin. The simulations
show that the different models and choice of initial conditions have small yet
subtle effects, with structures between simulations seeming to be rather
similar. To aid the eye we indicate the differences between matched particles
from $\Lambda$CDM to MG simulations with arrows/quivers. For the $\Lambda$CDM
ICs the MG simulations show more extreme features, _i.e._ denser clusters and
more emptier voids an effect driven by the larger clustering properties for
these models. However, for the MG ICs the differences are not as trivial, with
small displacements near clusters but large displacements between structures.
We interpret this as the constraints dictating the final location of the main
structures; leading to differences in the field between structures owing to
the model’s different formation histories.

To get a better comparison we compute the density and velocity field on a
$256^{3}$ grid using the triangular-shaped-cloud particle mass-assignment
scheme. We then isolate the cells within a $30\,h^{-1}{\rm Mpc}$ radius
[similar to 32] and compare the densities and radial velocities cell-by-cell
shown in Fig 6. Constrained simulation in $\Lambda$CDM, following the
procedure in this paper, have previously been shown to be consistent with
input peculiar velocities [see 32]. Therefore, by comparing to the
$\Lambda$CDM peculiar velocities and density fields we can test the accuracy
of the MG constrained simulations.

For MG constraints simulations produced with $\Lambda$CDM ICs we find the
density and radial velocities to be systematically shifted with low variance.
For the nDGP simulations produced with MG ICs we find consistent radial
velocities with $\Lambda$CDM. For $f(R)$, with MG ICs, we see a strong
negative correlation between radial velocities. Since the input are peculiar
velocities we expect these values to be consistent between models. This shows
the MG ICs for nDGP are performing as we expect and this method can be used
for future nDGP constrained simulations produced at higher resolution.
However, for $f(R)$ we find large discrepancies between the simulations
produced with MG ICs and find more consistent radial velocities with
$\Lambda$CDM ICs; although the latter comes at the cost of a systematic shift.
Since on large scales $f(R)$ has consistent growth functions as $\Lambda$CDM
this results shows that, in hindsight, our choice of scale was slightly
misguided and that we ought to have chosen a value at a larger scale, which
more closely reflect the scales that are more influential in structure
formation. For the time being, constrained simulations for $f(R)$ are better
produced with $\Lambda$CDM ICs rather than the approximate methods currently
used for MG ICs.

##  V Conclusions

In this paper, we extend the methodology for constructing constrained
simulations of the local Universe in $\Lambda$CDM to the MG models nDGP and
$f(R)$. We begin by describing the formal extensions to the methodology; the
input of different power spectra, correlation functions, and the numerical
calculation of growth functions. We use the ICeCoRe package [19], previously
used by CLUES [6] and HESTIA [8] for constructing the initial conditions for
constrained simulations in $\Lambda$CDM. Our implementation makes use of the
modular procedures to incorporate steps for MG. However, while this
implementation is completely compatible with the scale-independent growth
functions of nDGP, the implementation for $f(R)$ is incomplete. This is due to
the scale-dependence of the growth functions which are currently not
implemented in ICeCoRe. This issue arises due to a design choice in ICeCoRe to
constrain the density field from the displacement field and not the peculiar
velocities directly. The conversion from velocity to displacement includes a
division by the logarithmic growth function; for the scale-independent case
this is a simple conversion but for $f(R)$ this multiplication needs to be
carried out in Fourier space. For the peculiar velocity constraints applying
this conversion is impractical and since conducting this in a fully self-
consistent way would require rewriting significant portions of ICeCoRe we have
opted for an approximate scheme: taking an effective $f$ at a scale of
$k=2\pi/12$ $h{\rm Mpc}^{-1}$ (a real space scale of $R=12$ $h^{-1}{\rm Mpc}$)
roughly corresponding to the smoothing incurred from reconstructing the
density from peculiar velocities assuming linear theory.

We construct WF reconstructions of the density field from the CF2 peculiar
velocities for $\Lambda$CDM and MG nDGP and $f(R)$ models – the first such
construction for non-$\Lambda$CDM models. The WF reconstruction, on the whole,
show very similar features with subdominant variations owing to differences in
the MG linear power spectra; in nDGP features are more prominent at all scales
while for $f(R)$ this is limited to the largest amplitude features closest to
the origin (_i.e._ Local Group).

ICs for five constrained realisations were then generated with corresponding
ICs for random realisations with matching seeds. The WF reconstruction of the
mean density field show very similar reconstructions with subdominant
variations owing to differences in the MG linear power spectra. COLA MG
simulations were generated using both the $\Lambda$CDM and MG ICs.

Comparisons of the density field show that using $\Lambda$CDM ICs will
generate structures for MG that are more clustered and voids that are emptier,
this is incompatible with the constraints from peculiar velocity data and
shows the necessity for generating self-consistent MG ICs. Here we show that
the nDGP MG ICs reproduce consistent radial velocities with $\Lambda$CDM, a
property we expect since the constraints are from peculiar velocities. For
$f(R)$ the MG ICs are not fully self-consistent due to the scale-dependence of
the growth functions and the approximations used. As a result the $f(R)$
simulations with MG ICs produce radial velocities that are inconsistent with
$\Lambda$CDM. For $f(R)$ we find greater consistency with $\Lambda$CDM ICs. A
result showing that our choice of scale for the approximation of the growth
functions should have been larger, to reflect the scales which are more
important for structure formation. For the time being this means that
constrained simulations for $f(R)$ are more consistently produced with
$\Lambda$CDM ICs until a fully self-consistent scale-dependent methods can be
constructed.

In this paper we have extended the methodology for constructing constrained
simulations to MG. We have shown the importance of conducting this in a fully
self-consistent way (rather than using $\Lambda$CDM ICs), since MG imply
subtle but important differences in the reconstructed WF field and imply
different growth histories. Future work will look to extend this formalism to
fully incorporate the scale-dependence of models such as $f(R)$ without the
need to assume effective values for the growth rate. This will facilitate the
future production and study of high-resolution and hydrodynamic MG simulations
enabling the study of MG on properties of the local Universe on small scales.

###### Acknowledgements.

We thank Adi Nusser for providing useful comments and discussions. The
research conducted for this work was supported by the Polish National Science
Centre grants no: 2018/31/G/ST9/03388, 2018/30/E/ST9/00698,
2020/39/B/ST9/03494, 2020/38/E/ST9/00395, and by the Polish Ministry of
Science and Higher Education through grant DIR/WK/2018/12. NIL and SP
acknowledges support from the Deutsche Forschungs Gemeinschaft joint Polish-
German research project LI 2015/7-1. YH has been partially supported by the
Israel Science Foundation grant ISF 1358/18.

## References

  * Walker  _et al._ [1991] T. P. Walker, G. Steigman, D. N. Schramm, K. A. Olive, and H.-S. Kang, Primordial Nucleosynthesis Redux, ApJ 376, 51 (1991).
  * Planck Collaboration  _et al._ [2020] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J. P. Bernard, M. Bersanelli, _et al._ , Planck 2018 results. VI. Cosmological parameters, A&A 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Alam  _et al._ [2015] S. Alam, F. D. Albareti, C. Allende Prieto, F. Anders, S. F. Anderson, T. Anderton, B. H. Andrews, E. Armengaud, É. Aubourg, S. Bailey, S. Basu, J. E. Bautista, R. L. Beaton, T. C. Beers, C. F. Bender, _et al._ , The Eleventh and Twelfth Data Releases of the Sloan Digital Sky Survey: Final Data from SDSS-III, ApJS 219, 12 (2015), arXiv:1501.00963 [astro-ph.IM] .
  * Davis  _et al._ [1985] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White, The evolution of large-scale structure in a universe dominated by cold dark matter, ApJ 292, 371 (1985).
  * Hoffman and Ribak [1991] Y. Hoffman and E. Ribak, Constrained Realizations of Gaussian Fields: A Simple Algorithm, ApJ 380, L5 (1991).
  * Gottloeber  _et al._ [2010] S. Gottloeber, Y. Hoffman, and G. Yepes, Constrained Local UniversE Simulations (CLUES), arXiv e-prints , arXiv:1005.2687 (2010), arXiv:1005.2687 [astro-ph.CO] .
  * Sorce  _et al._ [2016] J. G. Sorce, S. Gottlöber, G. Yepes, Y. Hoffman, H. M. Courtois, M. Steinmetz, R. B. Tully, D. Pomarède, and E. Carlesi, Cosmicflows Constrained Local UniversE Simulations, MNRAS 455, 2078 (2016), arXiv:1510.04900 [astro-ph.CO] .
  * Libeskind  _et al._ [2020] N. I. Libeskind, E. Carlesi, R. J. J. Grand, A. Khalatyan, A. Knebe, R. Pakmor, S. Pilipenko, M. S. Pawlowski, M. Sparre, E. Tempel, P. Wang, H. M. Courtois, S. Gottlöber, Y. Hoffman, I. Minchev, _et al._ , The HESTIA project: simulations of the Local Group, MNRAS 498, 2968 (2020), arXiv:2008.04926 [astro-ph.GA] .
  * Jasche and Wandelt [2013] J. Jasche and B. D. Wandelt, Bayesian physical reconstruction of initial conditions from large-scale structure surveys, MNRAS 432, 894 (2013), arXiv:1203.3639 [astro-ph.CO] .
  * Wang  _et al._ [2016] H. Wang, H. J. Mo, X. Yang, Y. Zhang, J. Shi, Y. P. Jing, C. Liu, S. Li, X. Kang, and Y. Gao, ELUCID - Exploring the Local Universe with ReConstructed Initial Density Field III: Constrained Simulation in the SDSS Volume, ApJ 831, 164 (2016), arXiv:1608.01763 [astro-ph.CO] .
  * Lilow and Nusser [2021] R. Lilow and A. Nusser, Constrained realizations of 2MRS density and peculiar velocity fields: growth rate and local flow, MNRAS 507, 1557 (2021), arXiv:2102.07291 [astro-ph.CO] .
  * Zaroubi  _et al._ [1999] S. Zaroubi, Y. Hoffman, and A. Dekel, Wiener Reconstruction of Large-Scale Structure from Peculiar Velocities, ApJ 520, 413 (1999), arXiv:astro-ph/9810279 [astro-ph] .
  * Strauss and Willick [1995] M. A. Strauss and J. A. Willick, The density and peculiar velocity fields of nearby galaxies, Phys. Rep. 261, 271 (1995), arXiv:astro-ph/9502079 [astro-ph] .
  * Kourkchi  _et al._ [2020] E. Kourkchi, R. B. Tully, S. Eftekharzadeh, J. Llop, H. M. Courtois, D. Guinet, A. Dupuy, J. D. Neill, M. Seibert, M. Andrews, J. Chuang, A. Danesh, R. Gonzalez, A. Holthaus, A. Mokelke, D. Schoen, and C. Urasaki, Cosmicflows-4: The Catalog of $\sim$10,000 Tully-Fisher Distances, ApJ 902, 145 (2020), arXiv:2009.00733 [astro-ph.GA] .
  * Hu and Sawicki [2007] W. Hu and I. Sawicki, Models of f(R) cosmic acceleration that evade solar system tests, Phys. Rev. D 76, 064004 (2007), arXiv:0705.1158 [astro-ph] .
  * Dvali  _et al._ [2000] G. Dvali, G. Gabadadze, and M. Porrati, 4D gravity on a brane in 5D Minkowski space, Physics Letters B 485, 208 (2000), arXiv:hep-th/0005016 [hep-th] .
  * Doumler  _et al._ [2013a] T. Doumler, Y. Hoffman, H. Courtois, and S. Gottlöber, Reconstructing cosmological initial conditions from galaxy peculiar velocities - I. Reverse Zeldovich Approximation, MNRAS 430, 888 (2013a), arXiv:1212.2806 [astro-ph.CO] .
  * Doumler  _et al._ [2013b] T. Doumler, H. Courtois, S. Gottlöber, and Y. Hoffman, Reconstructing cosmological initial conditions from galaxy peculiar velocities - II. The effect of observational errors, MNRAS 430, 902 (2013b), arXiv:1212.2808 [astro-ph.CO] .
  * Doumler  _et al._ [2013c] T. Doumler, S. Gottlöber, Y. Hoffman, and H. Courtois, Reconstructing cosmological initial conditions from galaxy peculiar velocities - III. Constrained simulations, MNRAS 430, 912 (2013c), arXiv:1212.2810 [astro-ph.CO] .
  * Zaroubi  _et al._ [1995] S. Zaroubi, Y. Hoffman, K. B. Fisher, and O. Lahav, Wiener Reconstruction of the Large-Scale Structure, ApJ 449, 446 (1995), arXiv:astro-ph/9410080 [astro-ph] .
  * Tassev  _et al._ [2013] S. Tassev, M. Zaldarriaga, and D. J. Eisenstein, Solving large scale structure in ten easy steps with COLA, J. Cosmology Astropart. Phys 2013, 036 (2013), arXiv:1301.0322 [astro-ph.CO] .
  * Winther  _et al._ [2017] H. A. Winther, K. Koyama, M. Manera, B. S. Wright, and G.-B. Zhao, COLA with scale-dependent growth: applications to screened modified gravity models, J. Cosmology Astropart. Phys 2017, 006 (2017), arXiv:1703.00879 [astro-ph.CO] .
  * Vainshtein [1972] A. I. Vainshtein, To the problem of nonvanishing gravitation mass, Physics Letters B 39, 393 (1972).
  * Khoury and Weltman [2004] J. Khoury and A. Weltman, Chameleon cosmology, Phys. Rev. D 69, 044026 (2004), arXiv:astro-ph/0309411 [astro-ph] .
  * Gorski [1988] K. Gorski, On the Pattern of Perturbations of the Hubble Flow, ApJ 332, L7 (1988).
  * Nusser and Davis [1994] A. Nusser and M. Davis, On the Prediction of Velocity Fields from Redshift Space Galaxy Samples, ApJ 421, L1 (1994), arXiv:astro-ph/9309009 [astro-ph] .
  * Zel’Dovich [1970] Y. B. Zel’Dovich, Reprint of 1970A&A…..5…84Z. Gravitational instability: an approximate theory for large density perturbations., A&A 500, 13 (1970).
  * Tully  _et al._ [2013] R. B. Tully, H. M. Courtois, A. E. Dolphin, J. R. Fisher, P. Héraudeau, B. A. Jacobs, I. D. Karachentsev, D. Makarov, L. Makarova, S. Mitronova, L. Rizzi, E. J. Shaya, J. G. Sorce, and P.-F. Wu, Cosmicflows-2: The Data, AJ 146, 86 (2013), arXiv:1307.7213 [astro-ph.CO] .
  * Sorce  _et al._ [2017] J. G. Sorce, Y. Hoffman, and S. Gottlöber, Towards an optimal sampling of peculiar velocity surveys for Wiener Filter reconstructions, MNRAS 468, 1812 (2017), arXiv:1703.02045 [astro-ph.CO] .
  * Sorce [2015] J. G. Sorce, Minimization of biases in galaxy peculiar velocity catalogues, MNRAS 450, 2644 (2015), arXiv:1504.06968 [astro-ph.CO] .
  * Tully  _et al._ [2016] R. B. Tully, H. M. Courtois, and J. G. Sorce, Cosmicflows-3, AJ 152, 50 (2016), arXiv:1605.01765 [astro-ph.CO] .
  * Sorce  _et al._ [2014] J. G. Sorce, H. M. Courtois, S. Gottlöber, Y. Hoffman, and R. B. Tully, Simulations of the Local Universe constrained by observational peculiar velocities, MNRAS 437, 3586 (2014), arXiv:1311.2253 [astro-ph.CO] .