Is a sinusoidal signal the only one not distorted by an LTI system? In a signals and systems class in engineering school, I was shown a proof that a sinusoidal signal goes through a linear time-invariant (LTI) system undistorted; its shape and frequency are not changed - only its amplitude and phase (potentially) are. 

In other words, if the input is of the form $v_i = V\sin\left(\omega t + \phi\right)$, then in general the output is of the form $v_o = V^{\prime}\sin\left(\omega t + \phi^\prime\right)$.

This was an "aha! moment" for me in understanding the usefulness and power of the Fourier series/transform.

Now I wonder, however, is a sinusoidal signal the only type of (periodic) signal that goes through an LTI system undistorted? 


If so, I would love to be directed to a proof, and to hear any insight about why this is so.  
If not, is it possible to define a
new transform based on the alternative kernel?
Complex exponential function $x(t)=e^{i\omega_0 t}$ is the eigenfunction of an LTI system. That is, the output of an LTI system to this input is an scaled version of it. The scaling factor is the system function when evaluated at $\omega_0$. That is:
$$y(t)=H(\omega_0)e^{i\omega_0t} \tag{1}$$ 
The proof is not so difficult. Just substitute the input $x(t)=e^{i\omega_0 t}$ in the differential equation that describes the system, take the Fourier transform, and then simplify. Note that the given equation in the question is also derived from $(1)$.
The eigenvalues and eigenfunctions of LTI systems are given by the Laplace transform


In the continuous case :
not only the sines but also any exponential. With an input $f(t) = e^{s t}$ where $s  \in \mathbb{C}$ the ouput of your LTI system is 
$$f \ast h(t) = \int_{-\infty}^\infty e^{s (t-\tau)}h(\tau)d\tau = H(s) e^{st}$$
where $H(s)=\int_{-\infty}^\infty h(t) e^{-st}dt$ is the (bilateral) Laplace transform of the impulse response $h(t)$.
In the discrete case :
not only the sines but also any discrete exponential. With an input $f(n) = z^n$ where $z  \in \mathbb{C}$ the ouput is 
$$f \ast h(n) = \sum_{n=-\infty}^\infty z^{n-m}h(m) = H(z) z^n$$
where $H(z)=\sum_{-\infty}^\infty h(n)z^{-n}$ is the Z-transform of the impulse response $h(n)$.


And the only other eigenfunctions are the linear combination of such, if associated to the same eigenvalue $H(s),H(z)$.

Now in the context of signal processing, you'll accept slight modifications of the input signal (some sort of approximate eigenfunctions). Studying those is the purpose of filters theory, one of the most important problem being the approximation of the ideal low-pass filter.