  We present results of the first lattice QCD calculations of $B_c \to B_s$ and
$B_c \to B_d$ weak matrix elements. Form factors across the entire physical
$q^2$ range are then extracted and extrapolated to the physical-continuum limit
before combining with CKM matrix elements to predict the semileptonic decay
rates $\Gamma(B_c^+ \to B_s^0 \overline{\ell} \nu_{\ell}) = 26.2(1.2) \times
10^9 \,\text{s}^{-1}$ and $\Gamma(B_c^+ \to B^0 \overline{\ell} \nu_{\ell}) =
1.65(10) \times 10^9 \,\text{s}^{-1}$. The lattice QCD uncertainty is
comparable to the CKM uncertainty here. Results are derived from correlation
functions computed on MILC Collaboration gauge configurations with a range of
lattice spacings including 2+1+1 flavours of dynamical sea quarks in the Highly
Improved Staggered Quark (HISQ) formalism. HISQ is also used for the
propagators of the valence light, strange, and charm quarks. Two different
formalisms are employed for the bottom quark: non-relativistic QCD (NRQCD) and
heavy-HISQ. Checking agreement between these two approaches is an important
test of our strategies for heavy quarks on the lattice. From chained fits of
NRQCD and heavy-HISQ data, we obtain the differential decay rates $d\Gamma/ d
q^2$ as well as integrated values for comparison to future experimental
results.
