You reached a puzzling conclusion about $c_1(t) = c_2(t)$, and wonder whether you made a mistake in deriving them, or if the equality is indeed correct then how to explain it, perhaps by explicitly deriving one from the other.
I cannot tell whether it's possible to explicitly manipulate the double summation in $c_2(t)$ so as to convert it into the single summation of $c_1(t)$. It may be possible, I haven't tried. But sometimes it's not possible to explicitly show it, and instead, we have to rely on indirect evidence to confirm the equality and use the equality for our advantage if possible.
One such example is the Fourier transform of the constant $1$ :
$$ \mathcal{F}\{1\} = 2\pi \delta(\omega) \tag{0} $$
The equality in Eq.0 is not derived by explicitly evaluating the forward Fourier integral, but instead, deduced from the duality property of CTFT, and given that Fourier transform of $\delta(t)$ is $1$.
At the end of the analysis, we conclude that the validity of the equality $c_1(t) = c_2(t)$ is a consequence of the Nyquist sampling theorem;i.e., the truth of the equality is imposed by the sampling theorem, rather than a result of explicit algebraic manipulations of $c_2(t)$ into $c_1(t)$ or vice versa. And indeed, this is a useful side application of the theorem to prove that some equation is true which is very hard, if not impossible, to do so otherwise.
Let me show you, therefore, an indirect way of imposing the equality.
Let all signals $a(t),b(t)$, and $c(t)=a(t)b(t)$ are sufficiently bandlimited so that we can avoid aliasing.
Observe impulse train modulation relation:
$$ x_s(t) = x(t) \cdot \delta_T(t) ~ \cdot \cdot \cdot ~ \delta_T(t) = x(t) \cdot \delta_T(t) \tag{1}$$
where $ ~\delta_T(t) = \sum_n \delta(t-nT)$.
Also observe the interpolation relation :
$$ \left( x(t) \cdot \delta_T(t) \right) \star h(t) = x(t) \tag {2}$$
where $~h(t) = \text{sinc}(t/T) ~$ is the ideal lowpass brickwall interpolation filter.
We will use Eqs. 1 & 2 to derive alternative but equivalent expressions for the samples $c_s(t)$ of $c(t)$ and achieve the single and double sum versions $c_1(t)$ and $c_2(t)$.
The Double Sum :
$$
\begin{align}
c(t) &= c_s(t) \star h(t) \\
&= \left( c(t) \cdot \delta_T(t) \right) \star h(t) \\
&= \left( a(t) \cdot b(t) \cdot \delta_T(t) \right) \star h(t) &\text{Multiply then sample}\\
&= a(t)b(t) &\text{by Eq.2} \\
&= \left( \sum_n a[n] h(t-nT) \right) \left( \sum_m b[m] h(t-mT) \right) \\
&= \sum_n \sum_m a[n] b[n] ~ h(t-nT) ~ h(t-mT) \\
\end{align}
$$
The Single Sum :
$$
\begin{align}
c(t) &= c_s(t) \star h(t) \\
&= \left( a(t) \cdot b(t) \cdot \delta_T(t) \right) \star h(t) \\
&= \left( a(t) \cdot \delta_T(t) \cdot b(t) \cdot \delta(t) \right) \star h(t) &\text{by Eq.1} \\
&= \left( a_s(t) \cdot b_s(t) \right) \star h(t) &\text{Sample then multiply} \\
&= \left( \sum_n a[n] \delta(t-nT) \right) \left( \sum_m b[m] \delta(t-mT) \right) \star h(t) \\
&= \left( \sum_n \sum_m a[n]b[m] \delta((m-n)T) \cdot \delta(t-mT) \right) \star h(t) \\
&= \sum_n a[n] \left( \sum_m b[m] \left[ \delta((m-n)T)\delta(t-mT) \star h(t) \right] \right) \\
&= \sum_n a[n] \left( \sum_m b[m] \delta((m-n)T)) h(t-mT) \right) \\
&= \sum_n a[n]b[n] h(t-nT)
\end{align}
$$
The sampling relations are :
$$a_s(t) = a(t) \delta_T(t) = \sum_n a[n] \delta(t-nT) \tag{3}$$
$$b_s(t) = b(t) \delta_T(t) = \sum_m b[m] \delta(t-mT) \tag{4}$$
$$c_s(t) = c(t) \delta_T(t) = \sum_n c[n] \delta(t-nT) \tag{5}$$
We can get back $a(t),b(t)$,and $c(t)$ by ideal bandlimited interpolation of their samples $a[n],b[n]$,and $c[n]$ :
$$a(t) = a_s(t) \star h(t)= \sum_n a[n] \text{sinc}((t-nT)/T) \tag{6}$$
$$b(t) = b_s(t) \star h(t)= \sum_m b[m] \text{sinc}((t-mT)/T) \tag{7}$$
$$c(t) = c_s(t) \star h(t)= \sum_n c[n] \text{sinc}((t-nT)/T) \tag{8}$$
