If I have a continuous-time energy signal and multiply it with an impulse train, the frequency domain representation has infinite replicas. The total energy if we sum them up becomes infinite. What causes it to have infinite energy now?
I understand the derivation provided in many texts but I am trying to get intuitive understanding.
If you multiply a continuous-time finite energy signal $f(t)$ with an impulse train you get
$$\tilde{f}(t)=\sum_{n=-\infty}^{\infty}f(nT)\delta(t-nT)\tag{1}$$
where $T$ is the sampling interval and $\delta(t)$ is the Dirac delta impulse. Note that the "signal" $\tilde{f}(t)$ is a mathematical fiction, it cannot exist in practice, and it cannot even be evaluated for any value of $t$. Of course it is pretty safe to say that $\tilde{f}(t)$ is zero for $t\neq nT$, but you can't assign any value at $t=nT$. This comes from the fact that the Dirac delta impulse is no ordinary function but a generalized function or distribution, which only makes sense under an integral. So we can do calculations with $(1)$, we can take its Fourier transform, etc., but we cannot evaluate $(1)$ directly.
As for the energy of $\tilde{f}(t)$, since the square of $(1)$ is undefined (because the square of a distribution is undefined), we cannot compute the energy in the time domain. As you know, the spectrum of $(1)$ is the sum of shifted spectra of $f(t)$. Consequently, if you compute the energy in the frequency domain you get a diverging integral. So you could state that the energy of $(1)$ is definitely not finite. But the question remains if that has any meaning, given that $(1)$ is not even a function or signal in the ordinary sense.
Of course, the actual discrete-time signal $f[n]=f(nT)$, which is a sequence of finite numbers (the samples of $f(t)$), does have finite energy. If $f(t)$ is band-limited with maximum frequency $f_{max}<\frac{1}{2T}$ the energy of the sampled signal is equal to the energy of the continuous-time signal $f(t)$ (up to a constant):
$$E_f=\int_{-\infty}^{\infty}|f(t)|^2dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|^2d\omega\tag{2}$$
$$\begin{align}\sum_{n=-\infty}^{\infty}|f(nT)|^2&=\frac{T}{2\pi}\int_{-\pi/T}^{\pi/T}\left|\frac{1}{T}\sum_{k=-\infty}^{\infty}F\left(\omega-\frac{2\pi k}{T}\right)\right|^2d\omega\\&=\frac{1}{2\pi T}\int_{-\pi/T}^{\pi/T}|F(\omega)|^2d\omega=\frac{E_f}{T},\qquad\textrm{if } F(\omega)=0\textrm{ for }|\omega|>\pi/T\end{align}$$
A different way of showing this can be found in this answer.
When you have a more realistic sampling model which is an average of the signal over small time aperture, the overall infinite periodic spectra is low pass filtered, so the frequency images fall off as frequency increases.
Ideal delta function sampling is not physically realizable and as you add realistic physical effects the physical interpretations become more reasonable.
