So in the article about Green Functions in wikipedia definines
A Green's function, $G(x,s)$ of a linear differential operator $L=L(x)$ acting on distributions over a subset of the Euclidean space $\mathbb{R}^n$, at a point $s$, is any solution of $LG(x,s)=\delta(s-x)$ where $\delta$ is the Dirac delta function. This property of a Green's function can be exploited to solve differential equations of the form $Lu(x)=f(x)$
This definition seems to not require that the operator $L$ acts only on space, as in the same article, it later gives you a table of Green Functions for operators which depend on time, such as the heat operator.
The thing is that it is seems to have a hierarchy that is not explicitly formalized, at least not that I have understood. I would expect that the Green Function of an operator $L^\prime = \partial_s - L$ where $L$ is a elliptic operator is something of the sort $G(x,s;y,t)$ s.t
$$ (\partial_s - L)G(x,s;y,t) = \delta(t-s)\delta(x-y). $$
But in the case $L=\Delta$, the table gives us a
$$ G_t(x,y) = c e^{\frac{(x-y)^2}{4s}} $$ where $c$ is just a normalization constant. I understand that you can actually solve the PDE $$ \begin{cases} (\partial_s - \Delta)u(x,s) = f(x) , &(x,s)\in \mathbb{R} \times \mathbb{R}^+\\ u(x,0) = 0, & x \in \mathbb{R} \end{cases} $$ by using $u(x,s)= \int_{\mathbb{R}} G_s(x,y) f(y)dy$ and if I didn't do any wrong computation, we have that
$$ \int_{-\infty}^\infty (\partial_s - \partial_{xx})ce^{-\frac{(x-y)^2}{4s}} \phi(y) dy = c \int_{-\infty}^\infty \frac{e^{-\frac{(x-y)^2}{4s}}\phi(y) }{2s} dy, $$ which satisfies
$$ c\int_{-\infty}^\infty \frac{e^{-\frac{(x-y)^2}{4t}}\phi(y) }{2t} dy \longrightarrow f(x) \text{ as } t \to 0^+. $$
But I can't really see why one would have for $\phi(x,s)$ a test function $$ c \int_{0}^{\infty}\int_{-\infty}^\infty \frac{e^{-\frac{(x-y)^2}{4s}}\phi(y,s) }{2s} dy ds \stackrel{?}{=} \phi(x,0), $$ so does Green Function really behave like a Dirac in time?
So, how should I define the Green Function when there is time involved?
I face the same confusion when working with the eigenfunction expansion. For instance, for $L=\Delta$, it is easy to see that given a Hilbert basis of $L^2$ of eigenvectors $\{\phi_n\}_n$ with respective eigenvalues $\{\lambda_n\}_n$, then Green Function of $L$ is given by $$ G(x,y):= \sum_{n} \frac{\phi_n(x)\phi_n(y)}{\lambda_n}. $$ as for any test function $\phi$ $$ \int_D \phi(x) LG(x,y) dx = \sum_n \phi_n(y) \int \phi(x) \phi_n(x)dx = \sum_n \hat{\phi}_n \phi_n(y) = \phi(y). $$
But for the operator $L^\prime = \partial_s - \Delta$, Walsh's Book, affirm in page 323, that the Green Function of $L^\prime$ whith Neumann boundary conditions on $[0,L]$ is $$ G_s(x,y)= \sum_{k=1}^{\infty}\phi_n(x)\phi_n(y)e^{-\lambda_n s}, $$ where $\phi_n$ is still the eigenvector of $\Delta$ associated with $\lambda_n$. Again, I struggle to see how $G$ is being like a Dirac function in time. $$ \int_{0}^\infty \int_{0}^L\phi(y,s) L^\prime G_s(x,y)dxdt = -2 \int_{0}^\infty \int_0^L \sum_{k=1}^{\infty}\lambda_n\phi_n(x)\phi_n(y) \phi(y,s)e^{-\lambda_n s}dxdt \stackrel{?}{=} \phi(x,0). $$
So,
Is it matter of being in fact a slight different definition for Green Functions when the operator involves time? If so, what is the exact definition?
Or those Green functions actually behave like Dirac in time too? If so, why we only denote one parameter for time instead of the two parameter (as it is done for space)? Could anyone then give me a hint on how to complete the '$\stackrel{?}{=}$' steps?
First of all, some convenient notation. The points $x,y$ will live in space. The points $t,s$ in time.
For any possible operator $L_x$ (which differentiates w.r.t. $x$, not considering differences between time and space) a green function is defined as a $G(x,y)$ s.t. $$ L_x G(x, y) = \delta_{y}(x) \ \ $$ In this case $$ u_f(x) = \int dy G(x,y) f(y) $$ solves the inhomogeneous equation $$L_x u_f = f.$$ So far we do not care about initial conditions. On the other side, for an evolution equation: $$ (\partial_t - L_{x,t})u_f = f $$ you can use a more convenient notion of fundamental solution, namely one that solves for $t \le s $: $$ (\partial_t + L_{x,t})P(x,t; y,s) = 0 \\ P(x,s; y,s) = \delta_{y}(x) $$ as well as $$ (\partial_t - L^*_{y,s})P(x,t; y,s) = 0 \\ P(x,t; y,t) = \delta_{x}(y) $$ These equations have for example a nice probabilistic interpretation. In any case the fundamental solutions allow to solve the inhomogeneous problem This is linked to semigroup theory and ODEs: you just generalize Duhamel's formula). In the case of the heat equation the two equations coincide. Indeed the fundamental solution to the heat equation in symmetric (and the operator is self-adjoint).
A priori the two notions we introduced differ, but we can actually pass from one to the other. Take the fundamental solution to the heat equation $P_t(x).$ Then you the theory tells you that: $$ G(x,t;y,s) = \int_0^t dl \int dz P_{t-l}(x-z)\delta_y(x)\delta_s(l) = \begin{cases} 0 & t < s \\ P_{t-s}(x-y) & t \ge s \end{cases} $$
The latter is indeed a Green function for the heat operator. But as you can see it is more convenient to work with the heat kernel directly.