Link between LTI systems and square-integrable functions eigenfunctions I'm trying to study on my own a little bit of functional anaylsis. 
I try to summarize the main concepts (without being rigorous) that I learnt. Correct me where and if I'm wrong. 

In Hilbert spaces it is always possible to define a basis made of eigenfunctions and elements of the space can be rewritten using that basis. But how to find this basis? 
A self-adjoint operator $L$ on a space $H$, with the eigenvalue problem $Lx_{n}=\lambda x_{n}$, defines a basis of orthonormal eigenfunctions for the space objects and for the same operator expansion too. Moreover, there are particular linear self-adjoint operators called Sturm-Liouville operators that, defining an eigenvalue problem, generate too a basis for the corresponding space I'm working with. 

In particular, there are Sturm-Liouville "singular" problems that are defined over a space $H$ that is not "finite" but extends from $-\infty$ to $+\infty$: one of them, and for an engineer the most important, is $L^{2}(-\infty,+\infty)$. In this space, the operator $\frac{d^{2}}{dx^{2}}$ defines a Sturm Liouville singular problem and the solution of the eigenvalue problem is an infinite set of complex exponentials. From this theory we get the a continous spectrum for each signal of $L^{2}$ in terms of complex exponentials, i.e. a Fourier transform. I think that we can say that complex exponentials compose the spectrum of the operator $\frac{d^{2}}{dx^{2}}$ and of $L^{2}$ signals.

In engineering class, we met a particular class of operator called LTI (linear time-invariant) systems that have a peculiar characteristics: complex exponentials are eigenfunctions of these kind of operators.

My questions are: 

1) Is there some kind of link between the operator $\frac{d^{2}}{dx^{2}}$ spectrum and LTI operators? (PS: I'm not looking for the demonstration on wikipedia, but a justification or a reference coming from spectral theory)

2) I know that LTI systems can be characterized with the convolution integral: maybe the convolution integral defines the same eigenfunction basis with the corresponding eigenvalue problem?

3) I also know that some LTI systems can be written as $$\sum^{N}_{k=0}  a_{k}\frac{d^k}{dt^k}y(t)=\sum^{M}_{k=0} b_{k} \frac{d^k}{dt^k}x(t)$$, where basically we have TWO operators acting on two functions $x(t)$ and $y(t)$: these two functions can be written in the same eigenfunction basis? 

Sorry for the long post and thank you in advance if somebody can help me with these concepts or giving references! My usual references are Hanson-Yakovlev Operator Theory, Snider's PDE book and Oppenheim's Signals and Systems.
There are three kinds of Hilbert spaces: finite dimensional ones, separable infinite dimensional ones and inseparable infinite dimensional ones. The last category does not play an important role unless you practice quantum field theory, so you can normally ignore them. 

Two finite Hilbert spaces are equivalent if and only if they have the same dimension. Two separable Hilbert spaces are always equivalent. Separable Hilbert spaces have a countable basis.

Your question is about the spectral theorem. The spectral theorem comes in different flavours. Finite dimensional Hilbert spaces allow for a very simple formulation: Every self-adjoint operator $A$ has a real eigen-spectrum $\{ \lambda_k \in \mathbb{R} \vert k=1..K \}$ and decomposes into a finite sum of projectors $P_\lambda$ onto the eigen-subspace $E_\lambda$ associated with the eigenvalue $\lambda$.
$$A=\sum_{k=1}^K \lambda_k P_{\lambda_k}$$
The entire finite dimensional Hilbert decomposes into an orthogonal direct sum of the eigensubspaces.
$$\mathcal{H} = \bigoplus_{k=1}^K E_{\lambda_k} $$

This spectral decomposition only implies a unique orthogonal basis if all eigen-subspaces are one-dimensional in which case the corresponding operator is called non-degenerate. If you have a degenerate operator, you can choose an orthogonal basis for each eigen-subspace and get a orthogonal basis.

The more interesting case of a separable infinite dimensional Hilbert space makes things quite a bit more complicated. First, the concept of the spectrum of an operator needs to be generalised. In the finite dimensional case, the eigen-vector equation $A v = \lambda v$ for $v\neq0$ leads to the equivalent condition $\mathrm{ker}(A-\lambda)\neq\emptyset$. For infinite dimensional Hilbert spaces these two conditions are not equivalent and the requirement of a non-empty kernel is more general. Therefore a complex value $\lambda$ is defined to be part of the (generalised) spectrum of the operator $A$ exactly if $A-\lambda$ is singular. The associated kernel or nullspace is the generalised spectral subspace $E_\lambda:=\mathrm{ker}(A-\lambda)$.

The spectral theorem for separable infinite dimensional Hilbert spaces now states that for a compact self-adjoint operator $A$ we get a countable real spectrum and the spectral subspaces add up to the entire space. This is practically the situation we had before, only with infinite sums:
$$A=\sum_{k=1}^\infty \lambda_k P_{\lambda_k}$$
$$\mathcal{H}=\bigoplus_{k=1}^\infty E_{\lambda_k} $$
And again, an orthogonal basis of the Hilbert space is only implied if the operator is non-degenerate.

Compared to the finite dimensional case the big restriction is that the operator has to be compact. Unfortunately, many interesting operators are not compact, because they are not bounded. The spectral theory for non-compact operators is non-trivial and best avoided if possible. This is where Sturm-Liouville (SL) comes in. SL theory is a more specific version of the spectral theorem that works for a certain set of unbounded operators and also provides a method for calculating the spectral vectors. So if you can map an operator to a SL problem, you can make the same statements about the spectrum as for compact operators, even though the operator is unbounded.

Now, how does all this related to linear time-invariant (LTI) systems? Linear operators describe linear systems. Time-invariance can be added as an additional constraint. If you have a linear system operator $A$ and a unitary time-shift operator $S(\tau)$, then $A$ is time-invariant if and only if $AS(\tau)=S(\tau)A$, that is $A$ commutes with $S(\tau)$, for all real $\tau$. This is often written as $0=[A,S(\tau)]:=AS(\tau)-S(\tau)A$ with the commutator $[\cdot,\cdot]$. There always exists a common spectral basis for two commuting operators, that means both operators can be made diagonal in that basis. We will also need that an operator commutes with any function of that operator $[A,f(A)]=0$.

The time-shift operator can be expressed using a Taylor series generated by the exponential function $S(\tau)=\exp(\tau \frac{d}{dt})$. With what has been said above, it follows that LTI operators commute with $\frac{d}{dt}$ and therefore also with $\frac{d^2}{dt^2}$ and they share a spectral basis.

This should answer your first question. I will add more later to answer the two remaining ones.