How to realize Heaviside $\theta(t-t')$ and Dirac $\delta(t-t')$ as matrices in numerics? As is well known, single-particle Green's functions in the time domain might involve $$\theta(t-t')$$ for the retarded and advanced Green's functions. Sometimes, we also need $$\delta(t-t')$$ to express, e.g., instantaneous external potential. I am talking about the calculations based on various Dyson equations in the time domain, which often include time convolutions like
$$G_0\Sigma G \equiv \int d t_3 d t_4 G_0(t_1,t_3) \Sigma(t_3,t_4) G(t_4,t_2)$$
and even inversion operation like (convolution is understood)
$$[1 - G_0\Sigma]^{-1} G_0 .$$
One can numerically evaluate or solve them as discretized matrices. But how to correctly discretize $\theta(t-t')$ and $\delta(t-t')$ in a numerical realization? It's even not clear what to put for the '$1$' in the inversion above: the appropriately discretized $\delta(t-t')$ or just the identity matrix?
What is the good practice?
 A: I think in general it will depend on your application what the best solution is, since the answer will depend on the resolution of your simulation, how accurately you need to resolve local effects near the delta function, etc. A reasonable starting point would be to replace both functions with smooth approximations, such as
\begin{equation}
\delta(x-x') \rightarrow \frac{1}{\sqrt{2\pi \epsilon^2}} \exp\left(-\frac{(x-x')^2}{2\epsilon^2}\right)
\end{equation}
where $\epsilon$ is some small number, perhaps several grid points, and
\begin{equation}
\Theta(x-x') \rightarrow \frac{1}{2}\left(1+\tanh\left(\frac{x-x'}{\epsilon}\right)\right)
\end{equation}
where again, $\epsilon$ could be a few grid points.
As noted in the comments by Matteo (I will reproduce this here since comments can get deleted):

in the context of Green function often you have to deal with singularities of the form 1/(−) ( with ∈ℂ), where you often look for the imaginary part of this object. These singularities are often canceled by introducing a "Lorentzian broadening" Γ and modifying this function to 1/(−−Γ).

