  A study of $B_c^+\to J/\psi D_s^+$ and $B_c^+\to J/\psi D_s^{*+}$ decays
using 139 fb$^{-1}$ of integrated luminosity collected with the ATLAS detector
from $\sqrt{s} = 13$ TeV $pp$ collisions at the LHC is presented. The ratios of
the branching fractions of the two decays to the branching fraction of the
$B_c^+\to J/\psi \pi^+$ decay are measured: $\mathcal B(B_c^+\to J/\psi
D_s^+)/\mathcal B(B_c^+\to J/\psi \pi^+) = 2.76\pm 0.47$ and $\mathcal
B(B_c^+\to J/\psi D_s^{*+})/\mathcal B(B_c^+\to J/\psi \pi^+) = 5.33\pm 0.96$.
The ratio of the branching fractions of the two decays is found to be $\mathcal
B(B_c^+\to J/\psi D_s^{*+})/\mathcal B(B_c^+\to J/\psi D_s^+) = 1.93\pm0.26$.
For the $B_c^+\to J/\psi D_s^{*+}$ decay, the transverse polarization fraction,
$\Gamma_{\pm\pm}/\Gamma$, is measured to be $0.70\pm0.11$. The reported
uncertainties include both the statistical and systematic components added in
quadrature. The precision of the measurements exceeds that in all previous
studies of these decays. These results supersede those obtained in the earlier
ATLAS study of the same decays with $\sqrt{s} = 7$ and 8 TeV $pp$ collision
data. A comparison with available theoretical predictions for the measured
quantities is presented.
