# Implications for new physics from a novel puzzle in $\kern
1.79993pt\overline{\kern-1.79993ptB}{}{}^{0}_{(s)}\to
D^{(\ast)+}_{(s)}\\{\pi^{-},K^{-}\\}$ decays

Syuhei Iguro<EMAIL_ADDRESS>Department of Physics, Nagoya
University, Nagoya 464-8602, Japan  Teppei Kitahara
<EMAIL_ADDRESS>Institute for Advanced Research, Nagoya University,
Nagoya 464-8601, Japan  Kobayashi-Maskawa Institute for the Origin of
Particles and the Universe, Nagoya University, Nagoya 464-8602, Japan

###### Abstract

Recently, the standard model predictions for the $B$-meson hadronic decays,
$\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to D^{(\ast)+}K^{-}$ and
$\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}_{s}\to
D^{(\ast)+}_{s}\pi^{-}$, have been updated based on the QCD factorization
approach. This improvement sheds light on a novel puzzle in the $B$-meson
hadronic decays: there are mild but universal tensions between data and the
predicted branching ratios. Assuming the higher-order QCD corrections are not
huge enough to solve the tension, we examine several new physics
interpretations of this puzzle. We find that the tension can be partially
explained by a left-handed $W^{\prime}$ model, which can be compatible with
other flavor observables and collider bounds.

##  I Introduction

To test the standard model (SM) and search for physics beyond the SM,
precision measurements of meson decays, especially $B$-meson decays, have been
considerably investigated over the past 30 years. In the meantime, the
experimental uncertainty has been surprisingly reduced by experimentalists. On
the other hand, theorists have played an equally important role: several
approaches that can evaluate the QCD corrections have been invented, and the
SM predictions have been sharpened.

Very recently, SM predictions for several $B$-meson hadronic decays are
improved by Ref. Bordone:2020gao :

$\displaystyle\mathcal{B}(\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to
D^{+}K^{-})^{\rm exp}_{\rm SM}$
$\displaystyle=\left\\{\begin{array}[]{l}\left(1.86\pm 0.20\right)\times
10^{-4}\,,\\\ \left(3.26\pm 0.15\right)\times 10^{-4}\,,\end{array}\right.$
(3) $\displaystyle\mathcal{B}(\kern
1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to D^{\ast+}K^{-})^{\rm exp}_{\rm
SM}$ $\displaystyle=\left\\{\begin{array}[]{l}\left(2.12\pm 0.15\right)\times
10^{-4}\,,\\\ \left(3.27\,{}^{+0.39}_{-0.34}\right)\times
10^{-4}\,,\end{array}\right.$ (6) $\displaystyle\mathcal{B}(\kern
1.79993pt\overline{\kern-1.79993ptB}{}_{s}^{0}\to D^{+}_{s}\pi^{-})^{\rm
exp}_{\rm SM}$ $\displaystyle=\left\\{\begin{array}[]{l}\left(3.00\pm
0.23\right)\times 10^{-3}\,,\\\ \left(4.42\pm 0.21\right)\times
10^{-3}\,,\end{array}\right.$ (9) $\displaystyle\mathcal{B}(\kern
1.79993pt\overline{\kern-1.79993ptB}{}_{s}^{0}\to D^{\ast+}_{s}\pi^{-})^{\rm
exp}_{\rm SM}$ $\displaystyle=\left\\{\begin{array}[]{l}\left(2.0\pm
0.5\right)\times 10^{-3}\,,\\\ \left(4.3\,{}^{+0.9}_{-0.8}\right)\times
10^{-3}\,,\end{array}\right.$ (12)

where the upper numbers are the PDG averages of the experimental data PDG2020
, while the lower ones are the SM expectation values Bordone:2020gao . These
SM predictions are obtained by the QCD factorization (QCDF) Beneke:1999br ;
Beneke:2000ry ; Beneke:2001ev  at leading power in $\Lambda_{\rm QCD}/m_{b}$,
where the Wilson coefficients at next-to-next-to-leading logarithmic accuracy
are used Gorbahn:2004my . Compared to the previous estimations Huber:2016xod ,
the theoretical uncertainties are significantly reduced thanks to recent
developments in the $\kern 1.79993pt\overline{\kern-1.79993ptB}{}_{(s)}\to
D_{(s)}^{(\ast)}$ form factors including order $\mathcal{O}(1/m_{c}^{2})$
corrections within the framework of the heavy-quark expansion Jung:2018lfu ;
Bordone:2019vic ; Bordone:2019guc ; Iguro:2020cpg .

These hadronic channels are theoretically clean due to the absence of penguin
and annihilation topologies. Furthermore, resultant amplitudes are dominated
by the color-favored tree topology.

Above SM predictions deviate from the data at $5.6\sigma$ ($D^{+}K^{-}$),
$3.1\sigma$ ($D^{\ast+}K^{-}$), $4.6\sigma$ ($D^{+}_{s}\pi^{-}$), and
$2.4\sigma$ levels ($D^{\ast+}_{s}\pi^{-}$), respectively. Surprisingly, all
deviations are in the same direction and similar size. Note that
$\mathcal{B}(\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to
D^{+}\pi^{-})_{\rm SM}=(3.93{}^{+0.43}_{-0.42})\times 10^{-3}$ and
$\mathcal{B}(\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to
D^{\ast+}\pi^{-})_{\rm SM}=(3.45{}^{+0.53}_{-0.50})\times 10^{-3}$, which are
evaluated in Ref. Huber:2016xod , also deviate from the data,
$\mathcal{B}(\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to
D^{+}\pi^{-})^{\rm exp}=(2.52\pm 0.13)\times 10^{-3}$ and $\mathcal{B}(\kern
1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to D^{\ast+}\pi^{-})^{\rm
exp}=(2.74\pm 0.13)\times 10^{-3}$ PDG2020  at the $3.2\sigma$ and $1.4\sigma$
levels, respectively.

Within the SM, there are two possibilities that these tensions are alleviated.
The first possibility is an input value of $|V_{cb}|$. For $|V_{cb}|$, the
authors of Ref. Bordone:2020gao  use an average of the inclusive and exclusive
determinations in the $B$-meson semileptonic decays: $|V_{cb}|=(41.1\pm
0.5)\times 10^{-3}$ Bordone:2019vic ; Bordone:2019guc . If one adopts the
exclusive $|V_{cb}|$, $|V_{cb}|=(39.25\pm 0.56)\times 10^{-3}$ Amhis:2019ckw ,
amplitudes of the above processes are uniformly reduced by $4.5\%$. Note that
the exclusive $|V_{cb}|$, however, produces an additional $4.2\sigma$ level
tension in $\varepsilon_{K}$ Kim:2019vic . See also for a more recent
determination of the exclusive $|V_{cb}|$ using the full angular distribution
data Iguro:2020cpg .

Another possibility is higher-order QCD corrections. The next-to-leading power
and next-to-next-to-leading power corrections to the QCDF amplitudes are also
estimated by the same authors Bordone:2020gao , and the sizes of those
corrections to the amplitudes are evaluated as $\mathcal{O}(1)\%$.

The above puzzled situation could be resolved by introducing new physics
contributions to $b\to c\bar{u}q$ transitions, where $q=d$ and $s$.
Furthermore, it is shown that all ratios between these branching fractions are
consistent with data Bordone:2020gao . It clearly implies that the new physics
effects should be universal in $b\to c\bar{u}q$ transitions. Therefore, the
following questions are interesting: whether such a new physics is still
allowed by the other flavor constraints and by the hadron collider
constraints, and how much the tensions can be alleviated by a valid new
physics model. Below we will refer to this puzzle as $b\to c\bar{u}q$ anomaly.
In this Letter, we examine several new physics scenarios to explain the $b\to
c\bar{u}q$ anomaly.

##  II Framework

We consider the following effective Lagrangian to investigate new physics
contributions to $b\to c\bar{u}q$ processes:

$\displaystyle\mathcal{L}=-\frac{4G_{F}}{\sqrt{2}}\sum_{q}V_{cb}V_{uq}^{\ast}\sum_{i=1,2}C_{i}^{q}(\mu)\mathcal{Q}_{i}^{q}(\mu)\,,$
(13)

with the left-handed current-current operators in the CMM basis
Chetyrkin:1996vx ; Chetyrkin:1997gb ,

$\displaystyle\mathcal{Q}_{1}^{q}$
$\displaystyle=(\bar{c}_{L}\gamma^{\mu}T^{a}b_{L})(\bar{q}_{L}\gamma_{\mu}T^{a}u_{L})\,,$
(14) $\displaystyle\mathcal{Q}_{2}^{q}$
$\displaystyle=(\bar{c}_{L}\gamma^{\mu}b_{L})(\bar{q}_{L}\gamma_{\mu}u_{L})\,,$
(15)

where $q=d,s$. $T^{a}$ is the SU(3$)_{C}$ generator, and $V$ is the Cabibbo-
Kobayashi-Maskawa matrix Cabibbo:1963yz ; Kobayashi:1973fv . In our analysis,
we refrain from adding operators that are absent in the SM, e.g.,
$(\bar{c}_{L}b_{R})(\bar{q}_{L}u_{R})$. We will discuss this possibility in
the last section.

New physics contributions to the Wilson coefficients, $C_{1}^{q,{\rm NP}}$ and
$C_{2}^{q,{\rm NP}}$, become involved at the new physics scale $\Lambda$.
These values are modified by the renormalization-group (RG) evolution from
$\Lambda$ down to the hadronic scale $m_{b}$. The leading-order (LO) QCD RG
evolution is summarized in Appendix A. For instance, when $\Lambda=1\,{\rm
TeV}$, we obtain an evolution matrix as

$\displaystyle\begin{pmatrix}C_{1}^{\rm NP}(m_{b})\\\ C_{2}^{\rm
NP}(m_{b})\end{pmatrix}=\begin{pmatrix}1.36&-0.87\\\
-0.19&1.07\end{pmatrix}\begin{pmatrix}C_{1}^{\rm NP}(1\,{\rm TeV})\\\
C_{2}^{\rm NP}(1\,{\rm TeV})\end{pmatrix}\,.$ (16)

It is found that a universal destructive shift in the SM contributions is
favored in the $b\to c\bar{u}q$ anomaly Bordone:2020gao . The preferred size
is $\sim-17\%$, which corresponds to $C_{2}^{d,{\rm NP}}=C_{2}^{s,{\rm
NP}}=C_{2}^{\rm NP}$ and

$\displaystyle\frac{C_{2}^{\rm NP}(m_{b})}{C_{2}^{\rm SM}(m_{b})}=-0.17\pm
0.03\,.$ (17)

It is checked that such a new physics contribution is compatible with data of
the total decay rates of the $B$-mesons, $\tau_{B_{s}}/\tau_{B_{d}}$, and
$a_{d}^{fs}$ Bordone:2020gao ; Bobeth:2014rda ; Brod:2014bfa ; Lenz:2019lvd .
Another potentially strong constraint comes from the kaon hadronic decays
($s\to u\bar{u}d$). The ${\mathit{C}P}$-conserving parts of the isospin
amplitudes, $A_{I}=\langle(\pi\pi)_{I}|\mathcal{H}^{|\Delta S|=1}_{\rm
eff}|K\rangle$ for $I=0,2$, have been measured very precisely through all
$K\to\pi\pi$ data Kitahara:2016nld ; Blum:2015ywa

$\displaystyle\mbox{Re}A_{0}^{\rm exp}$ $\displaystyle=\left(3.3201\pm
0.0018\right)\times 10^{-7}\,{\rm GeV}\,,$ (18)
$\displaystyle\mbox{Re}A_{2}^{\rm exp}$ $\displaystyle=\left(1.4787\pm
0.0031\right)\times 10^{-8}\,{\rm GeV}\,.$ (19)

On the other hand, these theoretical predictions are

$\displaystyle\mbox{Re}A_{0}^{\rm SM}$ $\displaystyle=\left(2.99\pm
0.67\right)\times 10^{-7}\,{\rm GeV}\,,$ (20)
$\displaystyle\mbox{Re}A_{2}^{\rm SM}$ $\displaystyle=\left(1.50\pm
0.15\right)\times 10^{-8}\,{\rm GeV}\,,$ (21)

where the hadronic matrix elements are calculated by the lattice QCD
simulations Blum:2011ng ; Blum:2015ywa ; Bai:2015nea ; Ishizuka:2018qbn ;
Abbott:2020hxn . Although the $A_{2}$ amplitude is more sensitive to new
physics than $A_{0}$, we find that a $\pm 20\%$ new physics contribution to
the $s\to u\bar{u}d$ amplitude could be compatible with the data.

##  III Minimal flavor violation

First, we study the most simple possibility for new physics scenario: minimal
flavor violation (MFV) hypothesis DAmbrosio:2002vsn ; Isidori:2012ts . The
detailed calculations for this section can be found in Appendix B.

We examine a dimension-six operator,
$\mathcal{L}=1/(2\Lambda^{2})(\bar{Q}_{L}\gamma^{\mu}Q_{L})^{2}$, whose flavor
off-diagonal components are controlled by the quark Yukawa. In the quark mass-
diagonal basis, this operator produces $C_{2}^{q,{\rm
MFV}}(\Lambda)\sim-1/(2\sqrt{2}G_{F}\Lambda^{2})$. Then, the $b\to c\bar{u}q$
anomaly in Eq. (17) suggests $\Lambda\lesssim 0.49\,{\rm TeV}$.

Among the various flavor and collider constraints, a nonresonant dijet angular
distribution search in the LHC gives the most stringent constraint on this
scenario. The result is reported by the ATLAS collaboration at
$\sqrt{s}=13\,{\rm TeV}$ with $\int dt\mathcal{L}=37$ fb-1 Aaboud:2017yvp . We
interpret the result and obtain a $95\%$ C.L. exclusion limit as
$\Lambda<3.7\,{\rm TeV}$, which excludes the suggested $\Lambda\sim 0.49\,{\rm
TeV}$. From this collider constraint, we obtain a bound

$\displaystyle\frac{C_{2}^{\rm MFV}(m_{b})}{C_{2}^{\rm
SM}(m_{b})}\gtrsim-0.002\,.$ (22)

Hence, this scenario never explains the $b\to c\bar{u}q$ anomaly.

##  IV SU(2)$\times$SU(2)$\times$U(1) model

Next, we consider a new physics model that can produce a more convoluted
flavor structure. An extended electroweak gauge group
SU(2$)_{1}\times$SU(2$)_{2}\times$U(1$)_{Y}$ with heavy vectorlike fermions
produces heavy gauge bosons, $W^{\prime\pm}$ and $Z^{\prime}$, interacting
with the left-handed SM fermions with a nontrivial flavor structure
Langacker:1988ur ; Chivukula:2003wj ; Chiang:2009kb ; Boucenna:2016wpr ;
Boucenna:2016qad . These flavor structures are controlled by the number of
generations of the vectorlike fermions ($n_{\rm VF}$) and mixings between the
SM fermions and vectorlike fermions.

The heavy gauge boson interactions are Boucenna:2016qad

$\displaystyle\mathcal{L}=$
$\displaystyle+\frac{g_{ij}}{2}Z^{\prime}_{\mu}\bar{d}_{L}^{i}\gamma^{\mu}d_{L}^{j}-\frac{\left(VgV^{\dagger}\right)_{ij}}{2}Z^{\prime}_{\mu}\bar{u}_{L}^{i}\gamma^{\mu}u_{L}^{j}$
$\displaystyle-\frac{\left(Vg\right)_{ij}}{\sqrt{2}}W^{\prime+}_{\mu}\bar{u}_{L}^{i}\gamma^{\mu}d_{L}^{j}+{\rm
H.c.}\,,$ (23)

where $u_{L},d_{L}$ are the mass eigenstates, and a coupling $g_{ij}$ is
defined in the $d_{L}$ basis. In the following, we will take
$M_{W^{\prime}}=M_{Z^{\prime}}=M_{V}$ for simplicity. By integrating out
$W^{\prime\pm}$, new physics contribution $C_{2}^{q,{W^{\prime}}}$ is obtained
as

$\displaystyle C_{2}^{q,W^{\prime}}(M_{V})=$
$\displaystyle\frac{1}{4\sqrt{2}G_{F}M_{V}^{2}}\frac{(Vg)_{23}(Vg)_{1q}^{\ast}}{V_{cb}V_{uq}^{\ast}}\,.$
(24)

In order to generate an uniform shift in both $b\to c\bar{u}d$ and $b\to
c\bar{u}s$, a SM-like flavor structure in $(Vg)_{1q}$ is required, and hence
$g_{11}$ should be nonzero. When only $g_{11}$ is a nonzero entry in $g_{ij}$,
a dangerous $\bar{c}uZ^{\prime}$ flavor-changing neutral current is generated
and it is severely constrained by the $D$-meson mixing as
$|g_{11}|/M_{V}<{\mathcal{O}}(10^{-2})$ (TeV)-1 Golowich:2007ka . To evade
this bound, we follow the U(2$)^{3}$ flavor symmetry Barbieri:1995uv ;
Barbieri:2011ci  and take $g_{11}=g_{22}$ in $g_{ij}$ in the following
analyses. Then the bound from the $D$-meson mixing is significantly relaxed as
$|g_{11}|/M_{V}\lesssim 16\left(\rm{TeV}\right)^{-1}.$

Another flavor constraint comes from the $K\to\pi\pi$ data. By permitting a
$\pm 20\%$ new physics contribution to the Wilson coefficient of
$(\bar{u}_{L}\gamma^{\mu}d_{L})(\bar{s}_{L}\gamma_{\mu}u_{L})$ [see Eq. (21)],
we obtain

$\displaystyle|g_{11}|/M_{V}\lesssim 3.6\left(\rm{TeV}\right)^{-1}\,.$ (25)

Note that many types of diagrams contribute to $K\to\pi\pi$ decays, and
nonperturbative QCD plays an essential role there. Therefore, this bound is a
just reference value.

In addition to $g_{11}$, another nonzero entry of $g_{33}$ or $g_{23}$ is
necessary to produce $C_{2}^{q,W^{\prime}}$. Therefore, we consider the
following flavor texture:

$\displaystyle g_{ij}=\begin{pmatrix}g_{11}&0&0\\\ 0&g_{11}&g_{23}\\\
0&g_{23}&g_{33}\end{pmatrix}\,,$ (26)

and will discuss several scenarios in detail. We assume $g_{ij}$ is real for
simplicity. Note that when $g_{11}$ is $\mathcal{O}(1)$, production cross
sections of the heavy gauge bosons become considerably large in the hadron
collider, and hence we will mostly discuss the LHC constraints in each
subsection. To evade surveying a dedicated collider constraint for low-mass
region where the constraint would be more stringent, the mass range
$M_{V}>1\,{\rm TeV}$ is considered in our analysis.

###  IV.1 Scenario 1: $g_{11}$ and $g_{33}$

In this subsection, we take $g_{23}=0$ and consider a scenario of
$g_{ij}=\textrm{diag}(g_{11},\,g_{11},\,g_{33})$. Such a flavor structure can
be obtained from $n_{\rm VF}=1$. In this case, $(Vg)_{23}$ in Eq. (24) comes
from $V_{cb}g_{33}$. Since one has a factor of $V_{cb}$ just as the SM,
$\sqrt{|g_{11}g_{33}|}/M_{V}$ should be larger than $\mathcal{O}(1)\,{\rm
TeV}^{-1}$ to generate new physics contributions to $b\to c\bar{u}q$ processes
(see previous section). Furthermore, a relative sign between $g_{11}$ and
$g_{33}$ must be negative to produce the destructive interference with the SM
in the $b\to c\bar{u}q$ decays. A requirement of the $b\to c\bar{u}q$ anomaly
within 2 $\sigma$ level leads to

$\displaystyle 2.6\,({\rm
TeV})^{-1}\lesssim\sqrt{|g_{11}g_{33}|}/M_{V}\lesssim 3.8\,({\rm
TeV})^{-1}\,.$ (27)

Therefore, large couplings are necessary in this scenario.

First, let us examine the constraint from the $B_{s}$-meson mass difference
($\Delta M_{{s}}$). In this scenario, the dominant contribution comes from a
$W$–$W^{\prime}$ box diagram. We observed that the GIM mechanism still works
in this flavor structure, and obtain a simple formula for the $W$–$W^{\prime}$
box contribution to $\Delta M_{s}$,

$\displaystyle\frac{\Delta M_{{s}}^{W^{\prime}}}{\Delta M_{{s}}^{\rm
SM}}\simeq\eta^{\frac{2}{7}}\frac{2g_{11}g_{33}f^{\prime}(x_{t},x_{V})}{g_{W}^{2}f(x_{t})}\,,$
(28)

with $\eta=\alpha_{s}(M_{V})/\alpha_{s}(m_{W})$, $x_{t}=m_{t}^{2}/m_{W}^{2}$
and $x_{V}=M_{V}^{2}/m_{W}^{2}$, and $g_{W}$ is the weak coupling. The loop
functions are defined in Appendix C. We also have the same shift in
$B_{d}$-meson mixing, but it is less constrained because of its large
theoretical uncertainty. By imposing that the new physics contribution is
within $2\sigma$ uncertainty of $\Delta M_{{s}}^{\rm SM}$ Blanke:2018cya ;
DiLuzio:2019jyq , we obtain

$\displaystyle\sqrt{|g_{11}g_{33}|}/M_{V}\lesssim
1.7\left(\rm{TeV}\right)^{-1}\,.$ (29)

Although $\Delta M_{s}$ bound is incompatible with the $b\to c\bar{u}q$
anomaly in Eq. (27), we want to know how much this scenario can alleviate the
puzzle.

Figure 1:  Contours of $C_{2}^{\rm NP}(m_{b})/C_{2}^{\rm SM}(m_{b})$ are
presented. The puzzle can be explained at $2\sigma$ level in the yellow bands.
The blue and orange shaded regions are excluded by the dijet Aad:2019hjw ;
Sirunyan:2019vgj ; Sirunyan:2018xlo  and $t\bar{t}$ searches Sirunyan:2018ryr
; Aad:2020kop  at $95\%$ C.L., respectively. The regions above the dashed
lines are excluded by the single $t$ searches in the NWA (see text)
Sirunyan:2017vkm ; Aaboud:2018jux . Furthermore, the gray, red, green, and
purple shaded regions are constrained by $K\to\pi\pi$, $\Delta M_{s}$, $\Delta
M_{d}$, and $b\to s\gamma$, respectively. The dotted line indicates
$\Gamma_{V}/m_{V}$ and the red-hatched regions represent
$\Gamma_{V}/m_{V}>100\%$. Left: scenario 1. We take $g_{33}=-g_{11}$. Middle:
scenario 2. We take $g_{23}=-0.01(M_{V}/{\rm TeV})$. Right: scenario 3. We
take $M_{V}=1\,{\rm TeV}$ and $g_{11}=-3.6$.

Next, we consider constraints from resonant productions of the heavy gauge
bosons at the LHC. When $g_{11}$ and $g_{33}$ entries are nonzero,
$Z^{\prime}$ is produced via $pp\to q\bar{q}\to Z^{\prime}$ and also $pp\to
b\bar{b}\to Z^{\prime}$, while $W^{\prime\pm}$ is produced thorough $pp\to
q\bar{q}^{\prime}\to W^{\prime\pm}$ processes. When $M_{V}\gg m_{t}$, the
decay width of those particle is approximately given as,

$\displaystyle\Gamma_{V=W^{\prime},Z^{\prime}}\simeq\frac{2|g_{11}|^{2}+|g_{33}|^{2}}{16\pi}m_{V}\,.$
(30)

We find that relevant collider bounds come from dijet and $t\bar{t}$ searches.
The former provides the relevant bound for $|g_{11}|\gg|g_{33}|$, while the
latter for $|g_{11}|\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}|g_{33}|$.

Currently, both ATLAS and CMS collaborations reported upper limits on the
heavy dijet resonance cross section using the data of $\sim 140\,{\rm
fb}^{-1}$ Aad:2019hjw ; Sirunyan:2019vgj . Since $\mathcal{O}$(1) couplings
are necessary to relax the tension, the decay width can be not small.
Therefore, we adopt width-dependent limits on the cross section times the
dijet branching ratio. The broader the width is, the weaker the limits become
because a characteristic resonance peak is diluted. The search is robust up to
$\Gamma_{V}/M_{V}=20\%$ for $1.8$–$2.1\,{\rm TeV}$, and up to
$\Gamma_{V}/m_{V}=55\%$ for the heavier region Sirunyan:2019vgj . For the mass
range of $1$–$1.8\,{\rm TeV}$, we use an upper limit in Ref. Sirunyan:2018xlo
, where the narrow width approximation (NWA) is used. As for the heavy
$t\bar{t}$ resonance search, CMS reported the width-dependent limit using the
data of 36 ${\rm fb}^{-1}$ up to $\Gamma_{V}/M_{V}=30\%$ Sirunyan:2018ryr ,
while ATLAS reported the result using the data of 139 ${\rm fb}^{-1}$ in the
NWA Aad:2020kop .

We obtained the production cross section of $Z^{\prime}$ and $W^{\prime\pm}$
by rescaling the result in Refs. Sirunyan:2019vgj ; Hayreter:2019dzc , where
$\sigma(pp\to q\bar{q}^{\prime}\to W^{\prime+})+\sigma(pp\to
q\bar{q}^{\prime}\to W^{\prime-})\simeq 2\sigma(pp\to q\bar{q}\to Z^{\prime})$
is used Abe:2015uaa . The excluded regions from the dijet and $t\bar{t}$
searches are shown as the blue and orange shaded regions in Fig. 1 (left),
respectively.

We also show constraints from the single $t$ searches by using the data of
$\sim 36\,$fb-1 of CMS Sirunyan:2017vkm  and ATLAS Aaboud:2018jux : the
regions above the dashed lines in Fig. 1 (left) are excluded. Note that both
analyses assume the narrow resonance, and no study exists for broad
resonances.

Taking a conservative position, regions above the plateaus of the shaded areas
can not be excluded, where the corresponding $\Gamma_{V}/M_{V}$ exceeds the
maximum width shown in each experimental result: $\Gamma_{V}/M_{V}>30\%$ in
the $t\bar{t}$ search, and $\Gamma_{V}/M_{V}>55\%$ for $2.1$–$5\,{\rm TeV}$
and $\Gamma_{V}/M_{V}>20\%$ for $1.8$–$2.1\,{\rm TeV}$ in the dijet search.
The horizontal blue dashed lines are extrapolations obtained by assuming the
analysis of Ref. Sirunyan:2018xlo  is valid up to $\Gamma_{V}/M_{V}=20\%$, and
should be taken with more care. We note that limits from the dijet angular
distribution data, which are not considered here, would also depend on the
width-mass ratio and only contact interaction models are investigated
ATLAS:2012pu ; Aaboud:2017yvp . Further dedicated analysis would be necessary
to exclude such a broad width region.

The red-hatched regions represent $\Gamma_{V}>m_{V}$, where a particle picture
is no longer valid and one could not discuss any conclusive prediction.

Note that both our study and above experimental analyses have considered only
the $s$-channel productions of $W^{\prime}$ and $Z^{\prime}$, although there
are several $t$-channel contributions. Since the $t$-channel process does not
show a resonant nature, and there is a huge QCD $t$-channel background in the
dijet production, we suppose that inclusion of the $t$-channel processes in
the signal could not amplify the signal-to-noise ratio in the resonance
searches. Such $t$-channel contributions, which are insensitive to the width,
would be potentially accessible in the angular distribution search.

As long as we allow the broad width scenario, we find that the bound from
$\Delta M_{s}$ in Eq. (29) determines the maximal deviation of
$C_{2}^{W^{\prime}}/C_{2}^{\rm SM}$, which is independent of the ratio of
$g_{11}$ and $g_{33}$. For these reasons, we conclude
$C_{2}^{W^{\prime}}/C_{2}^{\rm SM}\gtrsim-0.05$ when $g_{23}=0$.

###  IV.2 Scenario 2: $g_{11}$ and $g_{23}$

For the second scenario, we set $g_{33}=0$ and consider $g_{11}$ and $g_{23}$
in Eq. (26). Such a flavor structure can be obtained when $n_{\rm VF}=2$. In
this scenario, the $b\to c\bar{u}q$ anomaly requires

$\displaystyle 0.54\,({\rm
TeV})^{-1}\lesssim\sqrt{|g_{11}g_{23}|}/M_{V}\lesssim 0.78\,({\rm
TeV})^{-1}\,.$ (31)

Although the size of the required coupling product is much smaller than the
previous scenario, a severe bound on $g_{23}$ comes from $\Delta M_{s}$, where
there is a tree-level $Z^{\prime}$ exchange diagram. We obtain

$\displaystyle\frac{\Delta M_{{s}}^{Z^{\prime}}}{\Delta M_{{s}}^{\rm
SM}}\simeq\eta^{\frac{2}{7}}\frac{16\pi^{2}g_{23}^{2}}{g_{W}^{4}(V_{ts}V_{tb}^{\ast})^{2}x_{V}x_{t}f(x_{t})}\,,$
(32)

and find that $g_{23}$ always gives a positive shift in $\Delta M_{{s}}$. The
constraint from $\Delta M_{{s}}$ is Blanke:2018cya ; DiLuzio:2019jyq

$\displaystyle|g_{23}|/M_{V}\lesssim 0.01\,({\rm TeV})^{-1}\,.$ (33)

Therefore, $g_{11}\gtrsim 30\,(M_{V}/{\rm TeV})\gg 4\pi$ is required by Eqs.
(31) and (33), which implies that the $b\to c\bar{u}q$ anomaly can not be
explained by this scenario.

In this scenario, $|g_{23}|\ll|g_{11}|$ should be satisfied. This simplifies
the collider constraints because the production cross section is controlled
only by $|g_{11}|$, and the heavy gauge bosons decay into jets with
$\mathcal{B}\simeq 1$. The constraints are shown in Fig. 1 (middle). We find
$C_{2}^{W^{\prime}}/C_{2}^{\rm SM}\gtrsim-0.01$, where
$g_{23}=-0.01(M_{V}/{\rm TeV})$ is taken.

###  IV.3 Scenario 3: $g_{11}$, $g_{23}$ and $g_{33}$

To see maximum value of $|C_{2}^{W^{\prime}}/C_{2}^{\rm SM}|$ in this model,
we combine the first and second scenarios: all $g_{11}$, $g_{23}$, and
$g_{33}$ are non-zero entries. The point of this scenario is that the severe
bound from $\Delta M_{s}$ can be turned off by

$\displaystyle\frac{\Delta M_{{s}}^{W^{\prime}}}{\Delta M_{{s}}^{\rm
SM}}+\frac{\Delta M_{{s}}^{Z^{\prime}}}{\Delta M_{{s}}^{\rm SM}}\sim 0\,,$
(34)

where the $W^{\prime}$ contribution is destructive and the $Z^{\prime}$ one is
constructive in $\Delta M_{s}$ (see previous subsections). We find, however,
that even if the $\Delta M_{s}$ bound is turned off, $g_{11}g_{33}$ is still
constrained from the $\Delta M_{d}$ as

$\displaystyle\sqrt{|g_{11}g_{33}|}/M_{V}\lesssim
2.3\,\left(\rm{TeV}\right)^{-1}\,.$ (35)

This bound restricts the possible $W^{\prime}$ contribution to the $b\to
c\bar{u}q$ processes. Also, we have checked a constraint from $b\to s\gamma$
data. We conclude that the $b\to s\gamma$ bound is less sensitive than $\Delta
M_{d}$, see Appendix D.

Since $|g_{23}|\ll|g_{11}|,|g_{33}|$ still holds in this scenario, the
collider constraints are almost the same as the scenario 1. We focus on a
parameter region that the all LHC constrains are evaded by the broad width of
the heavy gauge bosons. In Fig. 1 (right), $C_{2}^{W^{\prime}}/C_{2}^{\rm SM}$
is shown on $g_{23}$–$g_{33}$ plane by fixing $M_{V}=1\,{\rm TeV}$ and
$g_{11}=-3.6$ corresponding to the maximum value allowed by the $K\to\pi\pi$
data in Eq. (25). Eventually, we obtain

$\displaystyle\frac{C_{2}^{W^{\prime}}(m_{b})}{C_{2}^{\rm
SM}(m_{b})}\gtrsim-0.10\,.$ (36)

##  V Discussion

Motivated by a recent improvement of the SM predictions on $\kern
1.79993pt\overline{\kern-1.79993ptB}{}^{0}\to D^{(\ast)+}K^{-}$ and $\kern
1.79993pt\overline{\kern-1.79993ptB}{}^{0}_{s}\to D^{(\ast)+}_{s}\pi^{-}$, we
investigated the size of possible several new physics contributions to these
processes. In spite of severe bounds from the other flavor observables and the
LHC searches, we conclude that a $-10\%$ shift in the $b\to c\bar{u}q$
amplitude is possible by the left-handed $W^{\prime}$ model. Such a new
physics contribution can reduce the tension in the $b\to c\bar{u}q$ processes.

Since $g_{22}=g_{11}$ is a necessary condition, this model also produces new
physics contributions to $b\to c\bar{c}s$ processes with the same size
Jager:2017gal ; Jager:2019bgk . Although they, e.g., $B^{+}\to J/\psi K^{+}$,
have been measured precisely, the SM predictions suffer from large
nonfactorizable corrections Li:2006vq ; Liu:2010zh ; Li:2020app . We,
therefore, expect that the $b\to c\bar{c}s$ processes are less sensitive than
$b\to c\bar{u}q$.

It is unclear whether the new physics scalar operator can explain the $b\to
c\bar{u}q$ anomaly, but it is an interesting direction to consider it. For
instance, within a general two Higgs doublet model, a charged Higgs
interaction is Iguro:2017ysu

$\displaystyle{\cal
L}=-H^{+}\bar{u}^{i}(V\rho_{d}P_{R}-\rho_{u}^{\dagger}VP_{L})_{ij}d^{j}+{\rm
H.c.}\,,$ (37)

where $(V\rho_{d})_{23}$ is stringently constrained by $\Delta M_{s}$ via a
heavy neutral Higgs exchange, while $(\rho_{u}^{\dagger}V)_{23}$ is less
constrained by the flavor and collider observables Iguro:2018qzf ; Hou:2018zmg
. Therefore, a potentially large contribution to the $b\to c\bar{u}q$
processes would be expected.

## Acknowledgements

The authors thank Gauthier Durieux, Motoi Endo, Satoshi Mishima, Yael Shadmi,
Yotam Soreq and Michihisa Takeuchi for valuable comments and discussion on the
analysis. The work of S.I. is supported by the Japan Society for the Promotion
of Science (JSPS) Research Fellowships for Young Scientists, No.19J10980. S.I.
would like to thank the warm hospitality at KEK where he stayed during the
work. The work of T.K. is supported in part by the JSPS Grant-in-Aid for
Early-Career Scientists, No.19K14706. The work is also supported by the JSPS
Core-to-Core Program (Grant No.JPJSCCA20200002).

##  Appendix A Renormalization-group evolution

The LO RG evolution in the effective Lagrangian in Eq. (5) is given as
Chetyrkin:1996vx

$\displaystyle\frac{d\vec{C}(\mu)}{d\ln\mu}=\frac{\alpha_{s}(\mu)}{4\pi}\begin{pmatrix}-4&12\\\
\frac{8}{3}&0\end{pmatrix}\vec{C}(\mu)\,.$ (38)

According to Ref. Buras:1991jm , we obtain an analytic solution of the LO RG
evolution as

$\displaystyle\begin{pmatrix}C_{1}^{\rm NP}(m_{W})\\\ C_{2}^{\rm
NP}(m_{W})\end{pmatrix}=\begin{pmatrix}\frac{1}{3}\eta^{\frac{2}{7}}+\frac{2}{3}\eta^{-\frac{4}{7}}&\eta^{\frac{2}{7}}-\eta^{-\frac{4}{7}}\\\
\frac{2}{9}\eta^{\frac{2}{7}}-\frac{2}{9}\eta^{-\frac{4}{7}}&\frac{2}{3}\eta^{\frac{2}{7}}+\frac{1}{3}\eta^{-\frac{4}{7}}\end{pmatrix}\begin{pmatrix}C_{1}^{\rm
NP}(\Lambda)\\\ C_{2}^{\rm NP}(\Lambda)\end{pmatrix}\,,$ (39)

with $\eta=\alpha_{s}(\Lambda)/\alpha_{s}(m_{W})$. At the weak scale, the SM
contributions enter as Gorbahn:2004my

$\displaystyle C_{1}^{q}(m_{W})=15\frac{\alpha_{s}(m_{W})}{4\pi}+C_{1}^{q,{\rm
NP}}(m_{W})\,,\qquad C_{2}^{q}(m_{W})=1+C_{2}^{q,{\rm NP}}(m_{W})\,,$ (40)

and their RG evaluation from the weak scale to the hadronic scale is

$\displaystyle\begin{pmatrix}C_{1}(m_{b})\\\
C_{2}(m_{b})\end{pmatrix}=\begin{pmatrix}\frac{1}{3}\bar{\eta}^{\frac{6}{23}}+\frac{2}{3}\bar{\eta}^{-\frac{12}{23}}&\bar{\eta}^{\frac{6}{23}}-\bar{\eta}^{-\frac{12}{23}}\\\
\frac{2}{9}\bar{\eta}^{\frac{6}{23}}-\frac{2}{9}\bar{\eta}^{-\frac{12}{23}}&\frac{2}{3}\bar{\eta}^{\frac{6}{23}}+\frac{1}{3}\bar{\eta}^{-\frac{12}{23}}\end{pmatrix}\begin{pmatrix}C_{1}(m_{W})\\\
C_{2}(m_{W})\end{pmatrix}\,,$ (41)

with $\bar{\eta}=\alpha_{s}(m_{W})/\alpha_{s}(m_{b})$.

##  Appendix B Minimal flavor violation

In this section, we give the detailed calculations for the MFV scenario. In
the MFV hypothesis, the
SU(3$)_{Q_{L}}\times$SU(3$)_{U_{R}}\times$SU(3$)_{D_{R}}$ flavor symmetry is
introduced and it is broken only by the Yukawa interactions DAmbrosio:2002vsn
; Isidori:2012ts . Under this hypothesis, the flavor structure is the same as
the SM one: the flavor-changing neutral currents are automatically suppressed.
For the $b\to c\bar{u}q$ anomaly, we consider the following dimension-six
operator,

$\displaystyle\mathcal{L}=\frac{1}{2\Lambda^{2}}\left\\{\bar{Q}_{L}^{i}\left[\delta_{ij}+a(Y^{u}Y^{u{\dagger}})_{i\neq
j}\right]\gamma^{\mu}Q_{L}^{j}\right\\}^{2}\,,$ (42)

with $Y^{u}=V^{{\dagger}}\textrm{diag}(y_{u},y_{c},y_{t})$, and $a$ is a
dimensionless coupling. In the quark mass-diagonal basis ($u_{L}^{\rm
diag}=Vu_{L},~{}d_{L}^{\rm diag}=d_{L}$), this operator produces

$\displaystyle\mathcal{L}\simeq\frac{1}{\Lambda^{2}}(V_{cb}+ay_{t}^{2}V_{ts}^{\ast})V_{uq}^{\ast}\left(\bar{c}_{L}\gamma^{\mu}b_{L}\right)\left(\bar{q}_{L}\gamma_{\mu}u_{L}\right)\,.$
(43)

So, we obtain

$\displaystyle C_{2}^{q,{\rm
MFV}}(\Lambda)=-\frac{1}{\Lambda^{2}}\frac{\sqrt{2}}{4G_{F}}\left(1+ay_{t}^{2}\frac{V_{ts}^{\ast}}{V_{cb}}\right)\,.$
(44)

From the operator in Eq. (42), we also obtain a constraint from the
$B_{s}$-meson mass difference ($\Delta M_{{s}}$) as (cf., Ref.
Silvestrini:2018dos ),

$\displaystyle\left|\Lambda/a\right|\gtrsim 7.9\,{\rm TeV}\,,$ (45)

where the LO RG effect is taken into account Bagger:1997gg ,

$\displaystyle C_{LL}(m_{W})=\eta^{\frac{2}{7}}C_{LL}(\Lambda)\,,$ (46)

with $\eta=\alpha_{s}(\Lambda)/\alpha_{s}(m_{W})$, and the latest SM
estimation of $\Delta M_{s}$ is adopted Blanke:2018cya ; DiLuzio:2019jyq . We
required the new physics contribution to $\Delta M_{s}$ does not change the SM
prediction at $2\sigma$ level.

On the other hand, from the $b\to c\bar{u}q$ anomaly in Eq. (9), we find

$\displaystyle\Lambda\sim\sqrt{1-a}\left(0.43{}^{+5}_{-3}\right){\rm TeV}\,.$
(47)

Therefore, we reach a requirement for the anomaly:

$\displaystyle\Lambda\lesssim 0.49\,{\rm TeV}\,\textrm{~{}and~{}}\,|a|\lesssim
0.06\,.$ (48)

However, such a contact interaction can be probed by a non-resonant dijet
angular distribution search in the LHC. The result is reported by the ATLAS
collaboration at $\sqrt{s}=13\,{\rm TeV}$ with $\int dt\mathcal{L}=37$ fb-1
Aaboud:2017yvp . We interpret the result in terms of the operator in Eq. (42),
and obtain a $95\%$ CL exclusion limit,

$\displaystyle\Lambda<3.7\,{\rm TeV}\,\textrm{~{}and~{}}\,4.9\,{\rm
TeV}<\Lambda<8.3\,{\rm TeV}\,.$ (49)

This bound is clearly incompatible with Eq. (48). From this constraint, we
obtain a bound

$\displaystyle\frac{C_{2}^{\rm MFV}(m_{b})}{C_{2}^{\rm
SM}(m_{b})}\gtrsim-0.002\,.$ (50)

Therefore, this new physics scenario never explains the $b\to c\bar{u}q$
anomaly.

##  Appendix C Loop functions

The loop functions $f(x)$ and $f^{\prime}(x,y)$ in Eq. (20) are defined by
Lee:1998qq  (cf. Choudhury:2004bh )

$\displaystyle f(x)=$
$\displaystyle\frac{4-11x+x^{2}}{4(1-x)^{2}}-\frac{3x^{2}{\rm{ln}}\,x}{2(1-x)^{3}}\,,$
(51) $\displaystyle f^{\prime}(x,y)=$
$\displaystyle\frac{1}{4y(x-y)^{2}(1-x)^{2}}\biggl{[}(1-x)(4x^{2}+4y^{2}+5x^{2}y-8xy-4xy^{2}-x^{3})$
$\displaystyle-3x^{2}(x-2y+xy){\rm{ln}}\left(\frac{y}{x}\right)-\frac{3x(x-y)^{2}}{y-1}{\rm{ln}}\,y\biggr{]}\,,$
(52)

where $\lim_{y\to 1}f^{\prime}(x,y)=f(x)$. We note that Ref. Lee:1998qq
contains a typo in its Eq. (22), where $-x^{2}$ in the last term of the first
line in the arXiv version must be replaced by $-x^{3}$.

The loop functions $f_{\gamma}(x)$ and $f_{g}(x)$ in Eq. (58) are defined by
Grinstein:1990tj

$\displaystyle f_{\gamma}(x)$
$\displaystyle=\frac{3x^{3}-2x^{2}}{4(x-1)^{4}}\ln
x+\frac{-8x^{3}-5x^{2}+7x}{24(x-1)^{3}}\,,$ (53) $\displaystyle f_{g}(x)$
$\displaystyle=\frac{-3x^{2}}{4(x-1)^{4}}\ln
x+\frac{-x^{3}+5x^{2}+2x}{8(x-1)^{3}}\,.$ (54)

##  Appendix D $W^{\prime}$ and $Z^{\prime}$ contributions to $b\to s\gamma$

The effective Lagrangian for the $b\to s\gamma$ process is

$\displaystyle\mathcal{L}=\frac{G_{F}}{\sqrt{2}}V_{ts}^{\ast}V_{tb}\left[\sum_{i=1}^{6}C_{i}(\mu)\mathcal{Q}_{i}(\mu)+C_{7\gamma}(\mu)Q_{7\gamma}(\mu)+C_{8g}(\mu)Q_{8g}(\mu)\right]\,,$
(55)

with

$\displaystyle
Q_{7\gamma}=\frac{e}{8\pi^{2}}m_{b}\bar{s}\sigma^{\mu\nu}(1+\gamma_{5})bF_{\mu\nu}\,,\qquad
Q_{8g}=\frac{g_{s}}{8\pi^{2}}m_{b}\bar{s}\sigma^{\mu\nu}T^{a}(1+\gamma_{5})bG^{a}_{\mu\nu}\,,$
(56)

and the operators $\mathcal{Q}_{1}$–$\mathcal{Q}_{6}$ are defined in Ref.
Buras:1998raa .

By integrating out the heavy gauge bosons, we obtain

$\displaystyle C_{2}(M_{V})$
$\displaystyle\simeq\frac{g_{11}g_{33}}{g_{W}^{2}}\frac{m_{W}^{2}}{M_{V}^{2}}\,,$
$\displaystyle C_{3}(M_{V})$
$\displaystyle\simeq-\frac{g_{11}g_{23}}{2g_{W}^{2}V_{ts}^{\ast}}\frac{m_{W}^{2}}{M_{V}^{2}}\,,$
(57) $\displaystyle C_{7}(M_{V})$
$\displaystyle=\frac{g_{11}g_{33}}{g_{W}^{2}}\frac{m_{W}^{2}}{M_{V}^{2}}f_{\gamma}\left(\frac{m_{t}^{2}}{M_{V}^{2}}\right)\,,\qquad$
$\displaystyle C_{8}(M_{V})$
$\displaystyle=\frac{g_{11}g_{33}}{g_{W}^{2}}\frac{m_{W}^{2}}{M_{V}^{2}}f_{g}\left(\frac{m_{t}^{2}}{M_{V}^{2}}\right)\,,$
(58)

and remaining coefficients are set to zero at $\mu=M_{V}$. To obtain new
physics contributions at the hadronic scale, we solved the corresponding RG
evolution down to $\mu=m_{b}$ numerically:

$\displaystyle\frac{d\vec{C}(\mu)}{d\ln\mu}=\frac{\alpha_{s}(\mu)}{4\pi}\left(\hat{\gamma}^{(0){\rm
eff}}\right)^{T}\vec{C}(\mu)\,,\qquad\vec{C}=\left(C_{1},\,C_{2},\,\cdots,\,C_{6},\,C_{7}^{\rm
eff},\,C_{8}^{\rm eff}\right)\,,$ (59)

where the anomalous dimension matrix $\hat{\gamma}^{(0){\rm eff}}$ is given in
Refs. Ciuchini:1993fk ; Buras:1998raa . The $C_{7}^{\rm eff},\,C_{8}^{\rm
eff}$ are the effective Wilson coefficients which are required to cancel a
regularization scheme dependence Buras:1993xp . In this model, $C_{7}^{\rm
eff}(M_{V})=C_{7}(M_{V})$ and $C_{8}^{\rm eff}(M_{V})=C_{8}(M_{V})$.

Using $C_{7}^{\rm eff}(m_{b})$, we obtain a constraint from $b\to s\gamma$
data, where we required the new physics contributions are within a $2\sigma$
uncertainty range Misiak:2015xwa ; Endo:2017ums . The bound is sensitive to
$g_{23}$ which comes from the $Z^{\prime}$ contribution to $C_{3}(M_{V})$. For
$g_{23}=0$, we obtain

$\displaystyle\sqrt{|g_{11}g_{33}|}/M_{V}\lesssim 2.5\,\left({\rm
TeV}\right)^{-1}\,.$ (60)

This bound is significantly alleviated for $g_{23}/g_{33}>0$ region, while it
becomes stronger for $g_{23}/g_{33}<0$ region.

## References

  * (1) M. Bordone, N. Gubernari, T. Huber, M. Jung, and D. van Dyk, A puzzle in $\bar{B}_{(s)}^{0}\to D_{(s)}^{+}\\{\pi^{-},K^{-}\\}$ decays and extraction of the $f_{s}/f_{d}$ fragmentation fraction. arXiv:2007.10338.
  * (2) Particle Data Group Collaboration, Review of Particle Physics, Prog. Theor. Exp. Phys. 2020 (2020) 083C01.
  * (3) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, QCD factorization for $B\to\pi\pi$ decays: Strong phases and CP violation in the heavy quark limit, Phys. Rev. Lett. 83 (1999) 1914–1917 [hep-ph/9905312].
  * (4) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, QCD factorization for exclusive, nonleptonic B meson decays: General arguments and the case of heavy light final states, Nucl. Phys. B 591 (2000) 313–418 [hep-ph/0006124].
  * (5) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, QCD factorization in $B\to\pi K,\pi\pi$ decays and extraction of Wolfenstein parameters, Nucl. Phys. B 606 (2001) 245–321 [hep-ph/0104110].
  * (6) M. Gorbahn and U. Haisch, Effective Hamiltonian for non-leptonic $|\Delta F|=1$ decays at NNLO in QCD, Nucl. Phys. B 713 (2005) 291–332 [hep-ph/0411071].
  * (7) T. Huber, S. Kränkl, and X.-Q. Li, Two-body non-leptonic heavy-to-heavy decays at NNLO in QCD factorization, JHEP 09 (2016) 112 [arXiv:1606.02888].
  * (8) M. Jung and D. M. Straub, Constraining new physics in $b\to c\ell\nu$ transitions, JHEP 01 (2019) 009 [arXiv:1801.01112].
  * (9) M. Bordone, M. Jung, and D. van Dyk, Theory determination of $\bar{B}\to D^{(*)}\ell^{-}\bar{\nu}$ form factors at $\mathcal{O}(1/m_{c}^{2})$, Eur. Phys. J. C 80 (2020) 74 [arXiv:1908.09398].
  * (10) M. Bordone, N. Gubernari, D. van Dyk, and M. Jung, Heavy-Quark expansion for ${{\bar{B}}_{s}\rightarrow D^{(*)}_{s}}$ form factors and unitarity bounds beyond the ${SU(3)_{F}}$ limit, Eur. Phys. J. C 80 (2020) 347 [arXiv:1912.09335].
  * (11) S. Iguro and R. Watanabe, Bayesian fit analysis to full distribution data of $\bar{B}\to D^{(*)}\ell\bar{\nu}$: $|V_{cb}|$ determination and New Physics constraints, JHEP 08 (2020) 006 [arXiv:2004.10208].
  * (12) HFLAV Collaboration, Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2018. arXiv:1909.12524. 2019 updated results and plots available at https://hflav.web.cern.ch/.
  * (13) LANL-SWME Collaboration, 2019 Update on $\varepsilon_{K}$ with lattice QCD inputs, PoS LATTICE2019 (2019) 029 [arXiv:1912.03024].
  * (14) K. G. Chetyrkin, M. Misiak, and M. Munz, Weak radiative B meson decay beyond leading logarithms, Phys. Lett. B 400 (1997) 206–219 [hep-ph/9612313]. [Erratum: Phys.Lett.B 425, 414 (1998)].
  * (15) K. G. Chetyrkin, M. Misiak, and M. Munz, $|\Delta F|=1$ nonleptonic effective Hamiltonian in a simpler scheme, Nucl. Phys. B 520 (1998) 279–297 [hep-ph/9711280].
  * (16) N. Cabibbo, Unitary Symmetry and Leptonic Decays, Phys. Rev. Lett. 10 (1963) 531–533.
  * (17) M. Kobayashi and T. Maskawa, CP Violation in the Renormalizable Theory of Weak Interaction, Prog. Theor. Phys. 49 (1973) 652–657.
  * (18) C. Bobeth, U. Haisch, A. Lenz, B. Pecjak, and G. Tetlalmatzi-Xolocotzi, On new physics in $\Delta\Gamma_{d}$, JHEP 06 (2014) 040 [arXiv:1404.2531].
  * (19) J. Brod, A. Lenz, G. Tetlalmatzi-Xolocotzi, and M. Wiebusch, New physics effects in tree-level decays and the precision in the determination of the quark mixing angle $\gamma$, Phys. Rev. D 92 (2015) 033002 [arXiv:1412.1446].
  * (20) A. Lenz and G. Tetlalmatzi-Xolocotzi, Model-independent bounds on new physics effects in non-leptonic tree-level decays of B-mesons, JHEP 07 (2020) 177 [arXiv:1912.07621].
  * (21) T. Kitahara, U. Nierste, and P. Tremper, Singularity-free next-to-leading order $\Delta$S = 1 renormalization group evolution and $\epsilon_{K}^{\prime}/\epsilon_{K}$ in the Standard Model and beyond, JHEP 12 (2016) 078 [arXiv:1607.06727].
  * (22) T. Blum et al., $K\rightarrow\pi\pi$ $\Delta I=3/2$ decay amplitude in the continuum limit, Phys. Rev. D 91 (2015) 074502 [arXiv:1502.00263].
  * (23) T. Blum et al., The $K\to(\pi\pi)_{I=2}$ Decay Amplitude from Lattice QCD, Phys. Rev. Lett. 108 (2012) 141601 [arXiv:1111.1699].
  * (24) RBC, UKQCD Collaboration, Standard Model Prediction for Direct CP Violation in K$\rightarrow$$\pi$$\pi$ Decay, Phys. Rev. Lett. 115 (2015) 212001 [arXiv:1505.07863].
  * (25) N. Ishizuka, K. I. Ishikawa, A. Ukawa, and T. Yoshié, Calculation of $K\to\pi\pi$ decay amplitudes with improved Wilson fermion action in non-zero momentum frame in lattice QCD, Phys. Rev. D 98 (2018) 114512 [arXiv:1809.03893].
  * (26) RBC, UKQCD Collaboration, Direct CP violation and the $\Delta I=1/2$ rule in $K\to\pi\pi$ decay from the standard model, Phys. Rev. D 102 (2020) 054509 [arXiv:2004.09440].
  * (27) G. D’Ambrosio, G. F. Giudice, G. Isidori, and A. Strumia, Minimal flavor violation: An Effective field theory approach, Nucl. Phys. B 645 (2002) 155–187 [hep-ph/0207036].
  * (28) G. Isidori and D. M. Straub, Minimal Flavour Violation and Beyond, Eur. Phys. J. C 72 (2012) 2103 [arXiv:1202.0464].
  * (29) ATLAS Collaboration, Search for new phenomena in dijet events using 37 fb-1 of $pp$ collision data collected at $\sqrt{s}=$13 TeV with the ATLAS detector, Phys. Rev. D 96 (2017) 052004 [arXiv:1703.09127].
  * (30) P. Langacker and D. London, Mixing Between Ordinary and Exotic Fermions, Phys. Rev. D 38 (1988) 886.
  * (31) R. S. Chivukula, H.-J. He, J. Howard, and E. H. Simmons, The Structure of electroweak corrections due to extended gauge symmetries, Phys. Rev. D 69 (2004) 015009 [hep-ph/0307209].
  * (32) C.-W. Chiang, N. G. Deshpande, X.-G. He, and J. Jiang, The Family SU(2$)_{l}\times$ SU(2$)_{h}\times$ U(1) Model, Phys. Rev. D 81 (2010) 015006 [arXiv:0911.1480].
  * (33) S. M. Boucenna, A. Celis, J. Fuentes-Martin, A. Vicente, and J. Virto, Non-abelian gauge extensions for B-decay anomalies, Phys. Lett. B 760 (2016) 214–219 [arXiv:1604.03088].
  * (34) S. M. Boucenna, A. Celis, J. Fuentes-Martin, A. Vicente, and J. Virto, Phenomenology of an $SU(2)\times SU(2)\times U(1)$ model with lepton-flavour non-universality, JHEP 12 (2016) 059 [arXiv:1608.01349].
  * (35) E. Golowich, J. Hewett, S. Pakvasa, and A. A. Petrov, Implications of $D^{0}$ \- $\bar{D}^{0}$ Mixing for New Physics, Phys. Rev. D 76 (2007) 095009 [arXiv:0705.3650].
  * (36) R. Barbieri, G. R. Dvali, and L. J. Hall, Predictions from a U(2) flavor symmetry in supersymmetric theories, Phys. Lett. B 377 (1996) 76–82 [hep-ph/9512388].
  * (37) R. Barbieri, G. Isidori, J. Jones-Perez, P. Lodone, and D. M. Straub, $U(2)$ and Minimal Flavour Violation in Supersymmetry, Eur. Phys. J. C 71 (2011) 1725 [arXiv:1105.2296].
  * (38) M. Blanke and A. J. Buras, Emerging $\Delta M_{d}$ -anomaly from tree-level determinations of $|V_{cb}|$ and the angle $\gamma$, Eur. Phys. J. C 79 (2019) 159 [arXiv:1812.06963].
  * (39) L. Di Luzio, M. Kirk, A. Lenz, and T. Rauh, $\Delta M_{s}$ theory precision confronts flavour anomalies, JHEP 12 (2019) 009 [arXiv:1909.11087].
  * (40) ATLAS Collaboration, Search for new resonances in mass distributions of jet pairs using 139 fb-1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, JHEP 03 (2020) 145 [arXiv:1910.08447].
  * (41) CMS Collaboration, Search for high mass dijet resonances with a new background prediction method in proton-proton collisions at $\sqrt{s}=$ 13 TeV, JHEP 05 (2020) 033 [arXiv:1911.03947].
  * (42) CMS Collaboration, Search for narrow and broad dijet resonances in proton-proton collisions at $\sqrt{s}=13$ TeV and constraints on dark matter mediators and other new particles, JHEP 08 (2018) 130 [arXiv:1806.00843].
  * (43) CMS Collaboration, Search for resonant $\mathrm{t}\overline{\mathrm{t}}$ production in proton-proton collisions at $\sqrt{s}=13$ TeV, JHEP 04 (2019) 031 [arXiv:1810.05905].
  * (44) ATLAS Collaboration, Search for $t\bar{t}$ resonances in fully hadronic final states in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector. arXiv:2005.05138.
  * (45) CMS Collaboration, Search for heavy resonances decaying to a top quark and a bottom quark in the lepton+jets final state in proton–proton collisions at 13 TeV, Phys. Lett. B 777 (2018) 39–63 [arXiv:1708.08539].
  * (46) ATLAS Collaboration, Search for vector-boson resonances decaying to a top quark and bottom quark in the lepton plus jets final state in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, Phys. Lett. B 788 (2019) 347–370 [arXiv:1807.10473].
  * (47) A. Hayreter, X.-G. He, and G. Valencia, LHC constraints on $W^{\prime},~{}Z^{\prime}$ that couple mainly to third generation fermions. arXiv:1912.06344.
  * (48) T. Abe, T. Kitahara, and M. M. Nojiri, Prospects for Spin-1 Resonance Search at 13 TeV LHC and the ATLAS Diboson Excess, JHEP 02 (2016) 084 [arXiv:1507.01681].
  * (49) ATLAS Collaboration, ATLAS search for new phenomena in dijet mass and angular distributions using $pp$ collisions at $\sqrt{s}=7$ TeV, JHEP 01 (2013) 029 [arXiv:1210.1718].
  * (50) S. Jäger, M. Kirk, A. Lenz, and K. Leslie, Charming new physics in rare B-decays and mixing? Phys. Rev. D 97 (2018) 015021 [arXiv:1701.09183].
  * (51) S. Jäger, M. Kirk, A. Lenz, and K. Leslie, Charming New $B$-Physics, JHEP 03 (2020) 122 [arXiv:1910.12924].
  * (52) H.-n. Li and S. Mishima, Penguin pollution in the $B^{0}\to J/\psi K_{S}$ decay, JHEP 03 (2007) 009 [hep-ph/0610120].
  * (53) X. Liu, Z.-Q. Zhang, and Z.-J. Xiao, $B\to(J/\Psi,\eta_{c})K$ decays in the perturbative QCD approach, Chin. Phys. C 34 (2010) 937–943 [arXiv:0901.0165].
  * (54) Y.-Q. Li, M.-K. Jia, and Z. Rui, Revisiting nonfactorizable contributions to factorization-forbidden decays of $B$ mesons to charmonium, Chin. Phys. C 44 (2020) 113104 [arXiv:2006.07625].
  * (55) S. Iguro and K. Tobe, $R(D^{(*)})$ in a general two Higgs doublet model, Nucl. Phys. B 925 (2017) 560–606 [arXiv:1708.06176].
  * (56) S. Iguro and Y. Omura, Status of the semileptonic $B$ decays and muon g-2 in general 2HDMs with right-handed neutrinos, JHEP 05 (2018) 173 [arXiv:1802.01732].
  * (57) W.-S. Hou, M. Kohda, and T. Modak, Constraining a Lighter Exotic Scalar via Same-sign Top, Phys. Lett. B 786 (2018) 212–216 [arXiv:1808.00333].
  * (58) A. J. Buras, M. Jamin, M. E. Lautenbacher, and P. H. Weisz, Effective Hamiltonians for $\Delta S=1$ and $\Delta B=1$ nonleptonic decays beyond the leading logarithmic approximation, Nucl. Phys. B 370 (1992) 69–104. [Addendum: Nucl.Phys.B 375, 501 (1992)].
  * (59) L. Silvestrini and M. Valli, Model-independent Bounds on the Standard Model Effective Theory from Flavour Physics, Phys. Lett. B 799 (2019) 135062 [arXiv:1812.10913].
  * (60) J. A. Bagger, K. T. Matchev, and R.-J. Zhang, QCD corrections to flavor changing neutral currents in the supersymmetric standard model, Phys. Lett. B 412 (1997) 77–85 [hep-ph/9707225].
  * (61) K. Y. Lee and J. C. Lee, $B$ \- $\bar{B}$ mixings and rare $B$ decays in the top flavor model, Phys. Rev. D 58 (1998) 115001 [hep-ph/9803248].
  * (62) S. R. Choudhury, N. Gaur, A. Goyal, and N. Mahajan, $B_{d}$ \- $\bar{B}_{d}$ mass difference in little Higgs model, Phys. Lett. B 601 (2004) 164–170 [hep-ph/0407050].
  * (63) B. Grinstein, R. P. Springer, and M. B. Wise, Strong Interaction Effects in Weak Radiative $\bar{B}$ Meson Decay, Nucl. Phys. B 339 (1990) 269–309.
  * (64) A. J. Buras in Les Houches Summer School in Theoretical Physics, Session 68: Probing the Standard Model of Particle Interactions, pp. 281–539.  1998\.  hep-ph/9806471.
  * (65) M. Ciuchini, E. Franco, L. Reina, and L. Silvestrini, Leading order QCD corrections to $b\to s\gamma$ and $b\to sg$ decays in three regularization schemes, Nucl. Phys. B 421 (1994) 41–64 [hep-ph/9311357].
  * (66) A. J. Buras, M. Misiak, M. Munz, and S. Pokorski, Theoretical uncertainties and phenomenological aspects of $B\to X_{s}\gamma$ decay, Nucl. Phys. B 424 (1994) 374–398 [hep-ph/9311345].
  * (67) M. Misiak et al., Updated NNLO QCD predictions for the weak radiative B-meson decays, Phys. Rev. Lett. 114 (2015) 221801 [arXiv:1503.01789].
  * (68) M. Endo, et al., Gluino-mediated electroweak penguin with flavor-violating trilinear couplings, JHEP 04 (2018) 019 [arXiv:1712.04959].


# Program-Based Strategy Induction for Reinforcement Learning

Carlos G. Correa1<EMAIL_ADDRESS>
Thomas L. Griffiths2,3<EMAIL_ADDRESS>
Nathaniel D. Daw1,2<EMAIL_ADDRESS>
1Princeton Neuroscience Institute, Princeton University, Princeton, New Jersey
2Department of Psychology, Princeton University, Princeton, New Jersey
3Department of Computer Science, Princeton University, Princeton, New Jersey

###### Abstract

Typical models of learning assume incremental estimation of continuously-
varying decision variables like expected rewards. However, this class of
models fails to capture more idiosyncratic, discrete heuristics and strategies
that people and animals appear to exhibit. Despite recent advances in strategy
discovery using tools like recurrent networks that generalize the classic
models, the resulting strategies are often onerous to interpret, making
connections to cognition difficult to establish. We use Bayesian program
induction to discover strategies implemented by programs, letting the
simplicity of strategies trade off against their effectiveness. Focusing on
bandit tasks, we find strategies that are difficult or unexpected with
classical incremental learning, like asymmetric learning from rewarded and
unrewarded trials, adaptive horizon-dependent random exploration, and discrete
state switching.

Keywords: program induction; reinforcement learning; heuristics; strategy
discovery;

Classic models of reinforcement learning (RL) often assume that humans and
other animals incrementally estimate expected rewards or other continuously-
varying decision variables, and make noisy decisions based on them (?, ?, ?).
However, there is increased recognition that the learning strategies actually
adopted by humans (and even animals) deviate qualitatively from these models
and instead reflect more discrete, idiosyncratic, and perhaps explicit
heuristics: elaborating simple state-switching strategies like win-stay, lose-
shift (?, ?, ?) and reflecting memories of individual trials rather than
summaries (?, ?, ?, ?, ?). Despite individual attempts to capture these
behaviors by elaborating classic models, there is no general formalism
spanning the space of such strategies or predicting what variants organisms
adopt in particular situations, and how.

Recent work on computational approaches to strategy discovery has started to
engage with these questions. For example, meta-RL has been used to train
recurrent neural networks (RNNs) to have recurrent dynamics that implement a
learning process (?, ?, ?) which approximate decision-making strategies that
are either task-optimized (?, ?) or supervised to distill organisms’ learning
behavior (?, ?, ?, ?). However, this approach is lacking in a number of
respects. First, neural networks require considerable effort to interpret
(making them awkward for theory discovery) and are data- and compute-intensive
to optimize (making them unlikely candidates for a process model of how
organisms discover strategies). Relatedly, it is not clear that they represent
the most appropriate strategy space over which to conduct meta-learning for
these purposes. While the divergence of task-optimized RNNs from observed
behavior might be addressed by introducing resource costs (?, ?, ?), capacity-
limited RNNs have low-dimensional continuous state spaces that do not capture
the types of memory-dependent dynamics discussed above.

Accordingly, we propose instead to use Bayesian program induction to identify
program-structured strategies for meta-learning in RL problems. This builds on
approaches used to study cognitive representations (?, ?, ?, ?) and discover
machine learning algorithms (?, ?) that have been adapted to sequential
decision making (?, ?, ?). Our Bayesian framework has a natural interpretation
as a resource-rational approach (?, ?) by letting the cognitive cost of a
program correspond to its description length, so that simpler strategies are
more probable under the prior but more effective strategies are more likely to
result in optimal behavior under the likelihood. Importantly, we can explore
the relationship between previously-reported behavior and the strategies
resulting from different trade-offs between effectiveness and program
simplicity. By focusing on program-structured strategies for RL, we explore a
broader space of strategies than has previously been considered. These
programs are interpretable and, since they are small and can be extended with
modularity and reuse (?, ?, ?), may also lead to a process-level model of
strategy induction.

Using this approach, we identify a range of strategies in various bandit
tasks, accounting for some previously observed behavioral phenomena that
deviate from classical incremental learning. We find instances where
strategies focus exclusively on reward instead of omission, resembling
behavioral signatures in rodents (?, ?) and asymmetric learning more broadly
(?, ?). We show how this approach supports adaptive, horizon-specific
randomness in exploration, as previously observed (?, ?). We also find
strategies that switch between discrete decision states focused on either
exploiting known options or exploring new options, as in previous studies (?,
?). Our framework provides an interpretable alternative to existing methods
for strategy discovery, which we hope can provide new insights in other
domains.

## Program-Based Strategy Induction

We formulate strategies as programs and use sampling-based Bayesian inference
to find programs that are both simple and effective. We first introduce a
formalism for tasks and strategies, then describe our framework for inference,
and close with implementation details.

### Tasks and Strategies

We study finite-horizon tasks where at every time step $t$ an agent takes an
action $a_{t}$, makes an observation $o_{t}$, and receives a reward $r_{t}$,
resulting over time in a history of agent-environment interactions
$h_{t}=(a_{1:t},o_{1:t},r_{1:t})$. Observations may only provide partial
information about the state of the task (?, ?), so future observations and
rewards depend on the full history of agent-environment interactions,
$p(o_{t+1},r_{t+1}\mid h_{t},a_{t+1})$. Agents are modeled as $p(a_{t+1}\mid
h_{t})$, producing a task-specific distribution of histories $p(h_{t+1}\mid
h_{t})=p(a_{t+1}\mid h_{t})p(o_{t+1},r_{t+1}\mid h_{t},a_{t+1})$.

Strategies require performing two computations at every time step: updating
internal representations to incorporate new observations and rewards, and
using that updated internal representation to act. Formally, strategies
generate an updated memory $m_{t+1}$, based on information from the previous
trial (memory $m_{t}$, action $a_{t}$, observation $o_{t}$, and reward
$r_{t}$), which we refer to as the state update function $f$,

$m_{t+1}=f(m_{t},a_{t},o_{t},r_{t}).$ (1)

Choice behavior is driven by the updated memory $m_{t+1}$ and, for
convenience, information from the previous trial, by way of a policy function
$g$ that returns unnormalized log probabilities for each action,

$\displaystyle\log p(a_{t+1}\mid h_{t})$ $\displaystyle\propto
g(m_{t+1},a_{t},o_{t},r_{t}).$ (2)

A completely specified strategy $\pi=(m_{1},q_{1},f,g)$ also requires an
initial memory $m_{1}$ and unnormalized log policy $\log p(a_{1})\propto
q_{1}$ for the first time step. At subsequent time steps, memory is first
updated with Eq. 1 and then action is taken with Eq. 2. Agents following a
strategy are evaluated based on the expected sum of rewards they accrue over
$T$ steps,

$V(\pi)=\mathop{\mathbb{E}}_{h_{T}}\sum_{t=1}^{T}r_{t}.$

### Inference

We use Bayesian program induction to search the space of strategies with a
generative prior over strategies and a likelihood based on the effectiveness
of the strategy.

Table 1:  Primitive operations for strategies.  Primitives |  Description
---|---
Arithmetic, Logic |
0, …, 49 |  Integers from 0 to 49 (inclusive)
+, * |  Addition, multiplication
-, 1/(x) |  Negation, multiplicative inverse
<, == |  Less than, equals
&&, ||, ! |  And, or, negation
if(c,x,y) |  Returns x if condition c is true, y otherwise
Vectors |
vec_full(x) |  A vector filled with the value $x$
vec_n(x1,…,xn) |  A vector where the first $n$ entries are supplied and others are 0, e.g., vec_2(x,y)=[x,y,0,0]
v[i] |  Returns $\texttt{i}^{th}$ entry of v
assign(v,i,x) |  Updated copy of v, with v[i]=x
add_assign(v,i,x) |  Updated copy of v, with v[i]=v[i]+x
Inputs |
prev_action |  Previous action, $a_{t}$
reward |  Previous reward, $r_{t}$
state |  Memory from previous trial $m_{t}$ for $f$ or current trial $m_{t+1}$ for $g$
Action probabilities |
logit(l) |  For two-action tasks, $l=\log\frac{p(a=0)}{p(a=1)}$
softmax(w,v) |  Uses unnormalized log probabilities in $v$, scaled by $w$
action(a) |  Takes action $a$
argmax(v) |  Takes action with earliest, maximum value in $v$

Each component of a strategy is a program, a composition of the primitive
operations defined in Table 1. The primitives correspond to simple arithmetic,
logic, and vectors, as well as primitives specific to decision-making, such as
previous actions, rewards, and action distributions. For simplicity, we let
agent states be vectors of length 4. The space of programs consists of all
valid combinations of primitive operations that return the appropriate type
(either a vector for $m_{1}$ and $f$, or action probabilities for $q_{1}$ and
$g$), which corresponds to a context-free grammar. The prior, $p(\pi)$,
assumes expansions in the grammar are sampled uniformly at random, except
terminals (like integers and inputs) are 8 times more likely than non-
terminals in order to avoid excessive nesting. The integers, from 1 onward,
are distributed according to a geometric distribution with probability $0.5$.
We represent actions with integers. Invalid programs (non-integral or out of
bounds actions/indices, computing $1/0$) were assigned a value of $-\infty$.
Programs that compute starting values ($m_{1}$, $q_{1}$) have a restricted set
of primitives, excluding conditionals, vector indexing, and other input
values.

Following the standard approach for planning by inference (?, ?, ?, ?), we
formulate our likelihood by introducing a Bernoulli-distributed random
variable $\Omega$ that indicates the optimality of a strategy, so that

$\log p(\Omega=1\mid\pi)\propto\beta V(\pi),$

with $\beta$ a positive weight on the value. Notably, optimal strategies
maximize this likelihood. Combining the likelihood and prior, the posterior
probability that programs are optimal is

$\log p(\pi\mid\Omega=1)\propto\beta V(\pi)+\log p(\pi).$ (3)

The value weight $\beta$ controls the relative contribution of simplicity and
task performance to the prior, so that simpler solutions are preferred when
$\beta\to 0$ and more effective solutions are preferred when $\beta\to\infty$.

For inference, we use Markov chain Monte Carlo (MCMC) with Metropolis-Hastings
acceptance to obtain samples from the posterior, as in prior work (?, ?, ?).
We run five chains for $5\times 10^{5}$ steps per weight $\beta$, which are
$[100,300,1000,3000,10000,30000]$ for the ? (?) task and
$[10,30,100,300,1000,3000]$ for others. Proposal distributions are subtree-
regeneration (?, ?), and resampling a primitive (retaining arguments),
swapping arguments, inserting a subtree between a parent and child, and
deleting a subtree while retaining an ancestor (?, ?). A subset of
${m_{1},q_{1},f,g}$ was proposed to, with each program independently sampled
for inclusion (probability $\frac{1}{10}$ for $m_{1}$ and $q_{1}$ and
$\frac{1}{5}$ for $f$ and $g$), avoiding empty subsets by always including a
program sampled proportionally to inclusion probabilities. Our implementation
is built on the open-source Fleet library used in previous studies (?, ?).

Value was estimated from $2\times 10^{3}$ rollouts for the ? (?) task and
$10^{4}$ rollouts in other cases, using expected immediate rewards as well as
fixed environmental and agent seeds per chain to decrease variance. To
facilitate comparison across chains, we compare normalized value so that
chance behavior corresponds to 0 and behavior guided by an oracle (i.e.
complete knowledge of the reward probabilities) corresponds to 1. Pareto
frontiers were computed by finding the maximal strategies in the space defined
by normalized value and the prior, using the the top 100 strategies from each
chain.

## Exploring Strategies for a Simple Bandit Task

(a)

(b)  WSLS
---
⬇ m_1 = vec_2(0,1) q_1 = action(0) f = state g = argmax(assign(
state,prev_action,reward))

(c)  accumulator
---
⬇ m_1 = vec_full(0) q_1 = action(1) f = add_assign( state,
if(==(reward,prev_action),1,0), 1 ) g = argmax(state)

(d)  deterministic partial
---
accumulator
⬇ m_1 = vec_full(0) q_1 = action(1) f = add_assign( state,prev_action,reward)
g = argmax(state)

(e)  stochastic partial
---
accumulator
⬇ m_1 = vec_1(0) q_1 = logit(0) f = add_assign( state,prev_action,reward) g =
softmax(1,state)

Figure 1:  Strategies for the two-armed bandit with stationary, Bernoulli
rewards. a) The Pareto frontier of deterministic strategies, which are maximal
points in the space defined by the normalized value and prior. Points
correspond to individual strategies. One stochastic strategy described in the
text is included (unfilled point) for comparison. Value is normalized so that
chance behavior has a value of 0 and behavior guided by an oracle has a value
of 1. Also shown is the Bayes-optimal solution to the task (dotted line).
Example strategies are shown in b-e), marked in a), and described in the text.

We first survey the strategies discovered for a two-armed bandit, with a
finite horizon of 20 trials and stationary, Bernoulli reward probabilities
drawn from a uniform distribution. The simplicity of this domain makes it
possible to compute the value of strategies exactly. It is also possible to
compare to the optimal solution, which has a simple analytic solution that can
be reasonably computed for this horizon.

We first examine strategies when value is exactly computed, restricting
attention to deterministic policies for computational reasons. We focus on the
Pareto frontier defined by simplicity and utility, that is, solutions that are
local maxima in the space defined by simplicity and utility (Fig. 1(a)).
Overall, the model exposes a progression of more elaborate and more effective
strategies parameterized by this trade-off.

One simple solution we identify is win-stay, lose-shift (WSLS; Fig. 1(b)),
where an action is repeated if the previous trial was rewarded, or an
alternative action is selected if the previous trial was not rewarded (?, ?).
This solution has a trivial state update $f$ that simply passes a starting
vector onward, with the strategy instead implemented entirely in the policy
function $g$. For the sake of clarity, we briefly walk through the evaluation
of this strategy. When prev_action=0 and reward=1, the result of assign is
$[1,1]$. Since argmax returns the earliest maximum value, the subsequent
action is $0$, so winning results in repeating the action $0$. On the other
hand, if reward=0, the result of assign is $[0,1]$, so argmax would result in
the action $1$, meaning that a loss resulted in a switch. The program behaves
similarly for the other action. While a more human-readable implementation
might use a conditional, this implementation’s careful use of indices avoids
the need for a conditional and happens to be more probable under the prior.

A much more effective strategy accumulates the evidence indicated by rewards,
so that wins increase the preference for an arm and losses decrease it. The
algorithm in Fig. 1(c) implements this by incrementing the state at a
conditional index, so that wins for an arm or losses for the other arm result
in an increment as indexed by the arm. Since actions are selected by using
argmax, the arm with the greatest accumulated evidence is selected at any
time.

Other discovered strategies worse than the accumulator corresponded to partial
accumulators that only incorporated part of the reward feedback in either an
arm-specific or reward-specific way (like in Fig. 1(d)). Strategies better
than the accumulator were typically minor variants of it, like giving greater
weight to outcomes for the arm that was not initially selected.

For validation, we examined the strategies discovered when value was estimated
and stochastic policies were permitted. We found similar policies as above,
and use of stochastic policies that were often effectively deterministic due
to large values of inverse temperatures.

One question is whether the resource-rational trade-offs implied by this model
rationalize idiosyncratic behaviors in biological learning. One notable random
policy is the stochastic partial accumulator in Fig. 1(e), which shows a
modest increase in value over WSLS. Interestingly, this strategy ignores
losses. This can result in poor performance with deterministic behavior, since
choice is fixed to the first rewarding option (Fig. 1(d), marked in Fig.
1(a)). However, partial accumulation can perform effectively with stochastic
choice. This example is broadly consistent with findings of asymmetric
learning from positive and negative prediction errors, such as exclusively
learning from reward by rodents in a reversal learning task (?, ?). It also
echoes findings of optimistic learning biases in humans (?, ?) that
theoretical accounts (in a related task with counterfactual feedback) have
shown to be normative in some settings (?, ?).

To examine this more closely, we searched exhaustively over a small space of
strategies to see whether exclusive accumulation of either rewarded or
unrewarded trials was more effective. The space of strategies was a
generalization of the partial accumulators in Fig. 1, testing starting
conditions, inverse temperature, and whether policies were deterministic or
not. When rewards were accumulated, the best strategy we found was Fig. 1(e).
When unrewarded trials were accumulated (replacing reward with
+(-(1),reward)), the best performing strategy was deterministic and
implemented WSLS. Thus, within this limited class of strategies, a bias
towards positive information is rational.

(a)

(b)  stochastic accumulator
---
⬇ m_1 = vec_1(0) q_1 = logit(0) f = add_assign( state,prev_action,reward) g =
softmax($\blacksquare$,state)

(c)

(d)

(e)

(f)

Figure 2:  Adaptive random exploration, using a stochastic accumulator. a) The
Pareto frontier for solutions in the long horizon condition. Many solutions
simply accumulate rewards, but increase the determinism of policies (less
probable under prior) to achieve greater value. b) The discovered accumulator.
The “$\blacksquare$” is a placeholder for the inverse temperature. c)
Replotted from Fig. 2b in ? (?). Across conditions, there is a horizon-
specific adjustment of decision noise. d) The optimal horizon-specific inverse
temperature for the stochastic accumulator leads to more random behavior for
the long horizon condition. Horizon-specific inverse temperature was selected
to maximize Eq. 3 with $\beta=300$. Rewards were scaled by $\frac{1}{100}$. e)
Replotted from Fig. 3a in ? (?). In the long horizon condition, there is
horizon-specific adjustment of decision noise. f) Memory magnitude in the
accumulator grows over time, resulting in less randomness for later trials in
the long horizon condition. Same inverse temperature as d).

## Adaptive Random Exploration

A key feature of optimal behavior in bandit tasks is adaptive use of
exploration—exploring is important early on, but over time behavior favors the
best choice given previous observations. ? (?) showed that humans explore
adaptively in a bandit task, with increased preference for unfamiliar and less
valuable arms for a longer horizon and decreased preference for both for a
shorter horizon. While horizon-dependent exploration of unfamiliar arms is
predicted by normative accounts, adaptive random exploration is hard to
explain since optimal behavior is deterministic.

The task is a two-armed bandit with normal rewards that starts with four
forced-choice trials, to control the information participants have about each
arm, then proceeds to either one or six free-choice trials. Forced-choice
trials ignored action probabilities, and an additional primitive was added so
that strategies could test whether the previous trial was a forced-choice
trial. Strategies were found for each horizon condition independently and, to
facilitate comparing conditions, value was the average received reward on
free-choice trials.

Applying our framework to search for stochastic strategies for this task, we
primarily identified an accumulator along the Pareto front (long horizon shown
in Fig. 2(a), short horizon was similar), with varied inverse temperatures for
the softmax (Fig. 2(b), increasing with the value weight because the prior
disfavors larger constants). This strategy performs best with a high inverse
temperature, resulting in near-deterministic behavior. In this section, we
will use this example program to account for the adaptive random exploration
found by ?, focusing on the equal-information conditions where it is most
directly exposed.

One pattern of adaptive random exploration was between different time
horizons: on the first free-choice trial, participants were more random in the
long horizon condition (Fig. 2(c)). Focusing on parametrically varying the
strategy in Fig. 2(b), we found horizon-specific inverse temperatures for the
softmax that maximized the posterior (Eq. 3). The inverse temperature was
larger for the short horizon, resulting in more deterministic behavior
compared to the long horizon (Fig. 2(d)), consistent with the behavioral
findings.

This horizon-dependent randomness results from an interplay between the
accumulator’s memory dynamics and the representation cost of the policy.
Accumulating rewards instead of weighting them means that the magnitude of
entries in the agent’s state will increase over time, leading to greater
determinism. So, for the long horizon condition, the inverse temperature
primarily serves to make early choices more deterministic. Since integers are
geometrically-distributed, larger values of the constant softmax temperature
are less probable. This penalizes the determinism of a policy, resembling
information-theoretic approaches to penalizing policy complexity (?, ?, ?).
Combined, this means that a strategy specific to the long horizon condition
can afford to be more random early on, avoiding the penalty associated with
being deterministic.

The second pattern of adaptive random exploration identified by ? (?) is
within the long horizon condition: participant choice was more random when
there were more remaining trials (Fig. 2(e)). While these results incorporated
free-choice trials, they attempted to avoid a confound between information and
reward by focusing on cases where each arm had been selected an equal number
of times. Using the best inverse temperature for the long horizon, we saw
similar results (Fig. 2(f)). As described above, the accumulator’s entries
grow over time, resulting in a natural shift toward more deterministic
behavior.

(a)

(b)  ⬇ m_1 = vec_full(0) q_1 = softmax(0,vec_full(0)) f = state g = softmax(
reward, assign(state,prev_action,7) )
---

(c)

(d)  ⬇ m_1 = vec_1(4) q_1 = softmax(0,vec_1(0)) f = vec_1( +(1,state[reward]))
g = softmax(6, assign( vec_full(state[reward]), prev_action, 4 ))
---

(e)

Figure 3:  Discovering strategies with discrete decision states for a bandit
task with non-stationary reward. a) The state-based choice model in ? (?),
featuring distinct states where behavior is either more exploratory or
exploitative. Adapted from Fig. 1c in ? (?). b) A WSLS strategy that randomly
samples from actions after a loss. c) A state machine for WSLS strategy.
States correspond to distinct action distributions. Edge color indicates
action, dashed edges correspond to losses, and edges with action probability
of less than 1% were excluded. The initial conditions of the strategy were
modified to simplify the state machine. d) A more complex strategy that
exploits after a single win, but requires consecutive losses to switch. e)
State machine for d). Wins always lead to a state at top (as indicated), so
wins are excluded elsewhere. States with 4+ consecutive losses (at bottom)
were collapsed because they have similar action probabilities and identical
transitions.

## Discrete Decision States

We next examine strategies that emerge for learning in a bandit task with non-
stationary rewards, which is difficult to solve with the accumulation models
presented in previous sections. ? (?) found that both choice behavior and
neural activity from monkeys in this setting were captured by a slightly
generalized WSLS model that involved switching between distinct states:
exploring by random sampling, and exploiting by selecting a preferred action
(Fig. 3(a)). Quickly switching between several distinctive regimes of behavior
is difficult to model using error-based learning.

We study the task reported by ?, which is a 3-armed bandit with Bernoulli
rewards. The probability of reward for each arm is in $(0.1,0.2,\ldots,0.9)$.
After every trial, the probability has a 10% chance of drifting up or down by
$0.1$, while staying in bounds. We set the horizon to 500 trials.

Our framework identifies strategies with similar motifs as the one proposed by
?, so we qualitatively examine two in this section. The simpler program in
Fig. 3(b) implements WSLS where losses lead to random sampling among options—a
loss means reward=0, so softmax returns a uniform distribution. We also show
the program’s dynamics as a state machine in Fig. 3(c), where each state
corresponds to a unique distribution over actions, and the arrows show how
wins and losses change the policy. In the state machine, the explore state
results in random choice. Wins lead to an action-specific exploit state, and
losses lead back to the explore state. This strategy precisely mirrors the
schematic proposed by ? However, while their detailed results are consistent
with state switching, they also show that both neural and behavioral measures
somewhat violate the Markov conditional independence properties implied by
their WSLS model, suggesting the actual state space is more complicated. We
explore one such strategy next.

The strategy in Fig. 3(d) (state machine in Fig. 3(e)) elaborates both the
exploration and exploitation phases. In particular, when exploiting (top) it
switches only after three successive losses; for exploration (bottom) it tries
random arms until a win, but avoids the most recent choice. The code
implements this by counting consecutive losses in the $0^{th}$ entry of
memory, resetting the counter to 1 following a win, and using the counter to
inform choice probabilities. Importantly, this strategy retains the discrete,
state-like nature of the proposal in ? Similar to the bandit algorithms
discussed above, this strategy also responds more readily to positive
feedback, which is also consistent with theoretical findings (in a related
task) that choice-confirmation bias is adaptive when reward is non-stationary
(?, ?).

## Discussion

We have proposed a framework for strategy induction, using inferential methods
to identify strategies that are both simple and effective. By searching for
program-structured strategies, our approach is able to find strategy
representations that are more interpretable than results from other
approaches. We use this framework to examine strategies for RL tasks that
trade off between simplicity and effectiveness. We were able to discover
strategies consistent with previously observed behavior, like asymmetric
sensitivity to prediction errors, adaptive random exploration, and use of
discrete decision states.

We were able to justify adaptive, horizon-dependent random exploration in our
framework, as previously observed (?, ?). Within a fixed class of strategies,
our framework can cast this adaptive behavior as the result of a resource-
rational trade-off. Other approaches could explain this by penalizing the
divergence of policies from a uniform distribution (?, ?, ?, ?). However, our
approach goes beyond simple information-theoretic penalties on action
distributions since it can identify structured, task-specific strategies and
penalize them based on their description length.

One limitation of our current approach is that the primary cost of a program
is its representational cost. This neglects any execution-related costs the
program might incur, like memory use. Penalizing trial-specific computation
costs could encourage strategies that conditionally perform fewer computations
in some trials, which, for example, could result in habits which can be
adaptive under a speed-accuracy trade-off (?, ?). Future research could also
adjust the primitive operations available to reflect cognitive resource
limitations. For example, many kinds of memory are known to be subject to
capacity limits or noise, like working memory (?, ?) and action values (?, ?).

Our approach is a computational-level account (?, ?) of strategy induction,
leaving open questions about the algorithms underlying strategy discovery.
Future work could examine the process-level bias induced by varying the number
of samples taken in MCMC or for value estimation. Another question is how
library learning, which adds common subprograms as primitives (?, ?), could
accelerate inference. Library learning offers a simple way to perform
continual learning, which, by contrast, is an open challenge for neural
networks (e.g., requiring task-specific network dynamics to be identified and
frozen to avoid catastrophic interference in research by ?, ?).

While we have focused on resource-rational strategies for learning in a
handful of bandit tasks, future research could apply this framework to
identify interpretable strategies for sequential decision making in settings
like planning, problem solving, and other learning tasks. Future work could
also apply this framework to discover strategies directly from behavioral
data, as others have done with RNNs (?, ?, ?). We hope that the continued
development of strategy discovery methods can accelerate our understanding of
cognition by identifying novel features of behavior (?, ?) and avoid
misinterpretation of behavioral signatures (?, ?, ?).

## Acknowledgments

We thank Evan M. Russek, Flora Bouchacourt, and Mark K. Ho for helpful
discussions about early versions of this project. This research was supported
by grants from the U.S. Army Research Office (ARO W911NF-16-1-0474, awarded to
NDD), National Institute of Mental Health (R01MH135587, awarded to NDD), and
the NOMIS foundation (awarded to TLG). The funders had no role in study
design, data collection and analysis, decision to publish, or preparation of
the manuscript.

## References

  * Akam, Costa, DayanAkam et al. Akam, T., Costa, R., Dayan, P.  (2015).  Simple plans or sophisticated habits? State, transition and learning interactions in the two-step task.  _PLOS Computational Biology_ , _11_(12), 1–25.
  * Binz SchulzBinz Schulz Binz, M., Schulz, E.  (2022).  Modeling human exploration through resource-rational reinforcement learning.  In _Advances in Neural Information Processing Systems_ (Vol. 35).
  * Bornstein, Khaw, Shohamy, DawBornstein et al. Bornstein, A. M., Khaw, M. W., Shohamy, D., Daw, N. D.  (2017).  Reminders of past choices bias decisions for reward in humans.  _Nature Communications_ , _8_(1), 15958.
  * Collins CockburnCollins Cockburn Collins, A. G. E., Cockburn, J.  (2020).  Beyond dichotomies in reinforcement learning.  _Nature Reviews Neuroscience_ , _21_(10), 576–586.
  * Collins FrankCollins Frank Collins, A. G. E., Frank, M. J.  (2012).  How much of reinforcement learning is working memory, not reinforcement learning? A behavioral, computational, and neurogenetic analysis.  _European Journal of Neuroscience_ , _35_(7), 1024–1035.
  * Dezfouli BalleineDezfouli Balleine Dezfouli, A., Balleine, B. W.  (2012).  Habits, action sequences and reinforcement learning.  _European Journal of Neuroscience_ , _35_(7), 1036–1051.
  * Dezfouli, Griffiths, Ramos, Dayan, BalleineDezfouli et al. Dezfouli, A., Griffiths, K., Ramos, F., Dayan, P., Balleine, B. W.  (2019).  Models that learn how humans learn: The case of decision-making and its disorders.  _PLOS Computational Biology_ , _15_(6), 1–33.
  * Duan et al.Duan et al. Duan, Y., Schulman, J., Chen, X., Bartlett, P. L., Sutskever, I., Abbeel, P.  (2016).  RL2: Fast reinforcement learning via slow reinforcement learning.  _arXiv preprint arXiv:1611.02779_.
  * Duncan ShohamyDuncan Shohamy Duncan, K. D., Shohamy, D.  (2016).  Memory states influence value-based decisions.  _Journal of Experimental Psychology: General_ , _145_(11), 1420–1426.
  * Duncker, Driscoll, Shenoy, Sahani, SussilloDuncker et al. Duncker, L., Driscoll, L., Shenoy, K. V., Sahani, M., Sussillo, D.  (2020).  Organizing recurrent network dynamics by task-computation to enable continual learning.  In _Advances in Neural Information Processing Systems_ (Vol. 33).
  * Ebitz, Albarran, MooreEbitz et al. Ebitz, R. B., Albarran, E., Moore, T.  (2018).  Exploration disrupts choice-predictive signals and alters dynamics in prefrontal cortex.  _Neuron_ , _97_(2), 450–461.e9.
  * Ellis, Albright, Solar-Lezama, Tenenbaum, O’DonnellEllis et al. Ellis, K., Albright, A., Solar-Lezama, A., Tenenbaum, J. B., O’Donnell, T. J.  (2022).  Synthesizing theories of human language with Bayesian program induction.  _Nature Communications_ , _13_(1), 5024.
  * Ellis et al.Ellis et al. Ellis, K., Wong, C., Nye, M., Sablé-Meyer, M., Morales, L., Hewitt, L., … Tenenbaum, J. B.  (2021).  Dreamcoder: Bootstrapping inductive program synthesis with wake-sleep library learning.  In _Proceedings of the 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation_ (pp. 835–850).
  * Findling, Skvortsova, Dromnelle, Palminteri, WyartFindling et al. Findling, C., Skvortsova, V., Dromnelle, R., Palminteri, S., Wyart, V.  (2019).  Computational noise in reward-guided learning drives behavioral variability in volatile environments.  _Nature Neuroscience_ , _22_(12), 2066–2077.
  * Goodman, Tenenbaum, Feldman, GriffithsGoodman et al. Goodman, N. D., Tenenbaum, J. B., Feldman, J., Griffiths, T. L.  (2008).  A rational analysis of rule-based concept learning.  _Cognitive Science_ , _32_(1), 108–154.
  * Iigaya, Fonseca, Murakami, Mainen, DayanIigaya et al. Iigaya, K., Fonseca, M. S., Murakami, M., Mainen, Z. F., Dayan, P.  (2018).  An effect of serotonergic stimulation on learning rates for rewards apparent after long intertrial intervals.  _Nature Communications_ , _9_(1), 2477.
  * Ji-An, Benna, MattarJi-An et al. Ji-An, L., Benna, M. K., Mattar, M. G.  (2023).  Automatic discovery of cognitive strategies with tiny recurrent neural networks.  _bioRxiv preprint 2023.04.12.536629_.
  * Kaelbling, Littman, CassandraKaelbling et al. Kaelbling, L. P., Littman, M. L., Cassandra, A. R.  (1998).  Planning and acting in partially observable stochastic domains.  _Artificial Intelligence_ , _101_(1), 99–134.
  * Lai GershmanLai Gershman Lai, L., Gershman, S. J.  (2021).  Policy compression: An information bottleneck in action selection.  In K. D. Federmeier (Ed.), _The psychology of learning and motivation_ (Vol. 74, pp. 195–232).  Academic Press.
  * Lau GlimcherLau Glimcher Lau, B., Glimcher, P. W.  (2005).  Dynamic response-by-response models of matching behavior in rhesus monkeys.  _Journal of the Experimental Analysis of Behavior_ , _84_(3), 555–579.
  * Lefebvre, Summerfield, BogaczLefebvre et al. Lefebvre, G., Summerfield, C., Bogacz, R.  (2022).  A normative account of confirmation bias during reinforcement learning.  _Neural Computation_ , _34_(2), 307–337.
  * LevineLevine Levine, S.  (2018).  Reinforcement learning and control as probabilistic inference: Tutorial and review.  _arXiv preprint arXiv:1805.00909_.
  * Lieder GriffithsLieder Griffiths Lieder, F., Griffiths, T. L.  (2020).  Resource-rational analysis: Understanding human cognition as the optimal use of limited computational resources.  _Behavioral and Brain Sciences_ , _43_ , e1.
  * MarrMarr Marr, D.  (1982).  _Vision_.  W. H. Freeman.
  * Miller, Eckstein, Botvinick, Kurth-NelsonMiller et al. Miller, K. J., Eckstein, M., Botvinick, M. M., Kurth-Nelson, Z.  (2023).  Cognitive model discovery via disentangled RNNs.  _bioRxiv preprint 2023.06.23.546250_.
  * Moskovitz, Miller, Sahani, BotvinickMoskovitz et al. Moskovitz, T., Miller, K., Sahani, M., Botvinick, M. M.  (2023).  A unified theory of dual-process control.  _arXiv preprint arXiv:2211.07036_.
  * Ortega et al.Ortega et al. Ortega, P. A., Wang, J. X., Rowland, M., Genewein, T., Kurth-Nelson, Z., Pascanu, R., … Legg, S.  (2019).  Meta-learning of sequential strategies.  _arXiv preprint arXiv:1905.03030_.
  * PalminteriPalminteri Palminteri, S.  (2023).  Choice-confirmation bias and gradual perseveration in human reinforcement learning.  _Behavioral Neuroscience_ , _137_(1), 78–88.
  * Parker et al.Parker et al. Parker, N. F., Cameron, C. M., Taliaferro, J. P., Lee, J., Choi, J. Y., Davidson, T. J., … Witten, I. B.  (2016).  Reward and choice encoding in terminals of midbrain dopamine neurons depends on striatal target.  _Nature Neuroscience_ , _19_(6), 845–854.
  * Piray DawPiray Daw Piray, P., Daw, N. D.  (2021).  Linear reinforcement learning in planning, grid fields, and cognitive control.  _Nature Communications_ , _12_(1), 4942.
  * Plonsky, Teodorescu, ErevPlonsky et al. Plonsky, O., Teodorescu, K., Erev, I.  (2015).  Reliance on small samples, the wavy recency effect, and similarity-based learning.  _Psychological Review_ , _122_(4), 621–647.
  * Real, Liang, So, LeReal et al. Real, E., Liang, C., So, D., Le, Q.  (2020).  AutoML-zero: Evolving machine learning algorithms from scratch.  In H. D. III A. Singh (Eds.), _Proceedings of the 37th International Conference on Machine Learning_ (Vol. 119, pp. 8007–8019).
  * Rescorla WagnerRescorla Wagner Rescorla, R. A., Wagner, A. R.  (1972).  A theory of Pavlovian conditioning: Variations in the effectiveness of reinforcement and nonreinforcement.  In A. H. Black W. F. Prokasy (Eds.), _Classical conditioning II: Current research and theory_ (pp. 64–99).  New York: Appleton-Century-Crofts.
  * RobbinsRobbins Robbins, H.  (1952).  Some aspects of the sequential design of experiments.  _Bulletin of the American Mathematical Society_ , _58_(5), 527–535.
  * Sutton BartoSutton Barto Sutton, R. S., Barto, A. G.  (2018).  _Reinforcement learning: An introduction_.  Cambridge, MA: MIT Press.
  * Toussaint StorkeyToussaint Storkey Toussaint, M., Storkey, A.  (2006).  Probabilistic inference for solving discrete and continuous state Markov decision processes.  In _Proceedings of the 23rd International Conference on Machine Learning_ (pp. 945–952).
  * Wang et al.Wang et al. Wang, J. X., Kurth-Nelson, Z., Kumaran, D., Tirumala, D., Soyer, H., Leibo, J. Z., … Botvinick, M.  (2018).  Prefrontal cortex as a meta-reinforcement learning system.  _Nature Neuroscience_ , _21_(6), 860–868.
  * Wilson, Geana, White, Ludvig, CohenWilson et al. Wilson, R. C., Geana, A., White, J. M., Ludvig, E. A., Cohen, J. D.  (2014).  Humans use directed and random exploration to solve the explore–exploit dilemma.  _Journal of Experimental Psychology: General_ , _143_(6), 2074–2081.
  * Wingate, Goodman, Roy, Kaelbling, TenenbaumWingate et al. Wingate, D., Goodman, N. D., Roy, D. M., Kaelbling, L. P., Tenenbaum, J. B.  (2011).  Bayesian policy search with policy priors.  In _Proceedings of the Twenty-Second International Joint Conference on Artificial Intelligence - Volume Two_ (pp. 1565–1570).
  * Yang PiantadosiYang Piantadosi Yang, Y., Piantadosi, S. T.  (2022).  One model for the learning of language.  _Proceedings of the National Academy of Sciences_ , _119_(5), e2021865119.