We observe the process $\Upsilon(2S)\to D_s^{(*)+} D_{sJ}^{-}$ and continuum production $e^+e^- \to D_s^{(*)+} D_{sJ}^- $ at $\sqrt{s} = 10.52$ GeV (and their charge conjugates) using the data samples collected by the Belle detector at KEKB, where $D_{sJ}^-$ is $D_{s1}(2536)^-$ or $D^{*}_{s2}(2573)^-$. Both $D_{sJ}^-$ states are identified through their decay into $\bar{K}\bar{D}^{(*)}$. We measure the products of branching fractions ${\cal B}(\Upsilon(2S) \to D_{s}^{(*)+} D_{sJ}^-) {\cal B}(D_{sJ}^-\to \bar{K} \bar{D}^{(*)})$ and the Born cross sections $\sigma^{\rm Born}(e^+e^- \to D_{s}^{(*)+} D_{sJ}^-) {\cal B}(D_{sJ}^-\to \bar{K} \bar{D}^{(*)})$, and then compare the ratios $R_1 \equiv {\cal B}(\Upsilon(2S)\to D_{s}^{(*)+} D_{sJ}^-)/{\cal B}(\Upsilon(2S)\to\mu^{+}\mu^-)$ for $\Upsilon(2S)$ decays and $R_2 \equiv \sigma^{\rm Born}(e^+e^-\to D_{s}^{(*)+}D_{sJ}^-)/\sigma^{\rm Born}(e^+e^-\to \mu^{+}\mu^-)$ for continuum production. We obtain $R_1/R_2 = 9.7\pm 2.3 \pm 1.1$, $6.8 \pm 2.1 \pm 0.8$, $10.2 \pm 3.3 \pm 2.5$, and $3.4 \pm 2.1 \pm 0.5$ for the $D_s^+ D_{s1}(2536)^-$, $D_s^{*+} D_{s1}(2536)^-$, $D_s^+ D_{s2}^{*}(2573)^{-}$, and $D_s^{*+} D_{s2}^{*}(2573)^{-}$ final states in the $D_{sJ}^-\to K^{-} \bar{D}^{(*)0}$ modes, respectively. Therefore, the strong decay is expected to dominate in the $\Upsilon(2S)\to D_{s}^{(*)+}D_{sJ}^-$ processes. We also measure the ratios of branching fractions ${\cal B}(D_{s1}(2536)^-\to K_S^0 D^{*}(2010)^{-})/{\cal B}(D_{s1}(2536)^-\to K^{-} D^{*}(2007)^0) = 0.48 \pm 0.07 \pm 0.02$ and ${\cal B}(D_{s2}^{*}(2573)^- \to K_S^0 D^-)/{\cal B}(D_{s2}^{*}(2573)^- \to K^{-}D^0) = 0.49 \pm 0.10 \pm 0.02$, which are consistent with isospin symmetry. The second ratio is the first measurement of this quantity. Here, the first uncertainties are statistical and the second are systematic.
