The decays $$B_c^+ \rightarrow J/\psi D_s^+$$ and $$B_c^+ \rightarrow J/\psi D_s^{*+}$$ are studied with the ATLAS detector at the LHC using a dataset corresponding to integrated luminosities of 4.9 and 20.6 fb $$^{-1}$$ of pp collisions collected at centre-of-mass energies $$\sqrt{s} = 7$$ TeV and 8 TeV, respectively. Signal candidates are identified through $$J/\psi \rightarrow \mu ^+\mu ^-$$ and $$D_s^{(*)+}\rightarrow \phi \pi ^+(\gamma /\pi ^0)$$ decays. With a two-dimensional likelihood fit involving the $$B_c^+$$ reconstructed invariant mass and an angle between the $$\mu ^+$$ and $$D_s^+$$ candidate momenta in the muon pair rest frame, the yields of $$B_c^+ \rightarrow J/\psi D_s^+$$ and $$B_c^+ \rightarrow J/\psi D_s^{*+}$$ , and the transverse polarisation fraction in $$B_c^+ \rightarrow J/\psi D_s^{*+}$$ decay are measured. The transverse polarisation fraction is determined to be $$\Gamma _{\pm \pm }(B_c^+ \rightarrow J/\psi D_s^{*+})/\Gamma (B_c^+ \rightarrow J/\psi D_s^{*+}) = 0.38 \pm 0.23 \pm 0.07$$ , and the derived ratio of the branching fractions of the two modes is $$\mathcal {B}_{B_c^+ \rightarrow J/\psi D_s^{*+}}/\mathcal {B}_{B_c^+ \rightarrow J/\psi D_s^+} = 2.8 \,^{+1.2}_{-0.8} \pm 0.3$$ , where the first error is statistical and the second is systematic. Finally, a sample of $$B_c^+\rightarrow J/\psi \pi ^+$$ decays is used to derive the ratios of branching fractions $$\mathcal {B}_{B_c^+ \rightarrow J/\psi D_s^+}/\mathcal {B}_{B_c^+ \rightarrow J/\psi \pi ^+} = 3.8 \pm 1.1 \pm 0.4 \pm 0.2$$ and $$\mathcal {B}_{B_c^+ \rightarrow J/\psi D_s^{*+}}/\mathcal {B}_{B_c^+ \rightarrow J/\psi \pi ^+} = 10.4 \pm 3.1 \pm 1.5 \pm 0.6$$ , where the third error corresponds to the uncertainty of the branching fraction of $$D_s^+\rightarrow \phi (K^+K^-)\pi ^+$$ decay. The available theoretical predictions are generally consistent with the measurement.