The Rendering Equation
First thing first, the rendering equation in its solid angle formulation in vacuum is:
$$L(x,\omega) = L_e(x,\omega) + \int_{S^2} f(x, \omega_i, \omega_o) L(r(x,\omega), -\omega)|\cos\theta_x| d\sigma(\omega).$$
You forgot the emission function $L_e$ in your question. The solution $L:M\times S^2 \to \mathbb{R}$ (where $S^2=\{p\in\mathbb{R}^3\,:\, \|p\|_2=1\}$ is just the 2-sphere) is fully determined by the collection $(M, n, f, L_e)$, where $M\subset \mathbb{R}^3$ is the set of points forming the scene surfaces, $n:M\to S^2$ is the function that assigns a normal to every point of $M$, $f:M\times S^2 \times S^2 \to \mathbb{R}_+$ is the BSDF (it represents the light scattering properties of the materials of your surfaces), and $L_e:M\times S^2 \to \mathbb{R}_+$ is the emitted radiance function (it is non-zero on your light sources). It is fairly natural that your geometry $(M,n$), the materials $f$, and the light sources $L_e$ will determine the steady state of the radiance $L$ in the scene.
The function $r:M\times S^2 \to M\cup (\infty\cdot S^2)$ is the ray-tracing function that returns the first intersection point $r(x,\omega)$ along the ray with origin $x$ and direction $\omega$, if there is no intersection you can formally claim that it returns the intersection with some sphere at infinity (e.g. your environment map). Typically the rendering equation is formulated in terms of the incident radiance $L_i(x, \omega)$, however, in vacuum we have $L_i(x,\omega) = (GL)(x,\omega) = L(r(x,\omega),-\omega)$ - the radiance that arrives at $x$ from direction $\omega$ is equal to the radiance leaving $r(x,\omega)$ in direction $-\omega$, which is also fairly intuitive. The term $|\cos\theta_x| = |n(x)\cdot \omega|$ is due to Lambert's cosine law, but it can be derived also by just integrating a differential $2$-density on the sphere (think surface flux integrals from calculus 3). The $d\sigma(\omega)$ is the area measure on the sphere, i.e. if you take the standard spherical coordinates parametrization $x = \sin\theta\cos\phi$, $y=\sin\theta\sin\phi$, $z=\cos\theta$ and your rotate them so that $z$ is aligned with $n(x)$, then you may rewrite $d\sigma(\omega) = |\sin\theta| d\theta d\phi$ and thus:
$$L(x,\omega) = L_e(x,\omega) + \int_{0}^{2\pi}\int_0^{\pi} f(x, \omega_i, \omega_o) L(r(x,\omega), -\omega)|\cos\theta| |\sin\theta|d\theta d\phi.$$
Formal Solution
For brevity one typically writes the above in operator form $L = L_e +TL$ where $T:(M\times S^2\to \mathbb{R})\to (M\times S^2\to \mathbb{R})$ is the operator that stands for the integral. If we define the scattering operator $K$ and the propagation operator $G$:
$$(Kg)(x,\omega) = \int_{S^2}f(x,\omega,\omega_o)g(x,\omega)d\sigma_{x}^{\perp}(\omega), \quad (Gh)(x,\omega) = h(r(x,\omega),-\omega),$$
then $T=K\circ G$.
Since we have the identity $L=L_e+TL$ we can just plug it in itself once yielding $L=L_e + T(L_e + TL) = L_e+TL_e+T^2L$. Repeating this process until infinity yields the series $L = \sum_{k=0}^{\infty}T^kL_e$. Formally $\sum_{k=0}^{\infty}T^k$ is the solution operator for the rendering equation:
\begin{align}
L = L_e +TL \implies (I-T)L = L_e \implies L = (I-T)^{-1}L_e \implies (I-T)^{-1} = \sum_{k=0}^{\infty}T^k.
\end{align}
The latter is the Neumann series for $T$. While $\|T\|<1$ is a sufficient condition for convergence (i.e. if the brdf is energy conserving), it is not a necessary condition - for specific choices of $L_e$ the above series may be convergent even if $(I-T)^{-1}$ is unbounded. Note that $f,L,L_e$ are actual functions (i.e. $T$ maps functions to functions). In practice one may want to model specular reflection, point and directional lights etc., which according to "many treatments" cannot be modeled with classical functions and require distribution. Veach discusses this in his thesis (e.g. in some appendices on specular reflection and refraction), and even has a nice table 8.3 showing some examples.
Veach's Extension to Distributions
If you go to 5.A.1, Veach has a short note that $\int_{\mathbb{R}} f(x)\delta(x-x_0)\,dx$ is just symbolic notation for the evaluation functional $\Lambda_{x_0}(f) = f(x_0)$. This symbolic expression could also be framed as $\int_{\mathbb{R}} f(x)\delta(x-x_0)dx := \int_{\mathbb{R}} f(x)\,d\delta_{x_0}(x)$ where $\delta_{x_0}$ is the Dirac measure. Note however that $\delta$ is not the Radon-Nikodym derivative of this measure, i.e. the expression $\int_{\mathbb{R}} f(x)\delta(x-x_0)dx$ is still purely symbolic if not abuse of notation. When you write $$\int_{S^2} f(x,\omega,\omega_o)L_i(x,\omega)\,d\sigma_x^{\perp}(\omega) = \int_{S^2} \delta_{\sigma_x^{\perp}}(\omega-\omega_r)\rho(x,\omega,\omega_o)L_i(x,\omega)\,d\sigma_x^{\perp}(\omega),$$
the second expression is purely symbolic and means
\begin{align}
\int_{S^2} \delta_{\sigma_x^{\perp}}(\omega-\omega_r)\rho(x,\omega,\omega_o)L_i(x,\omega)\,d\sigma_x^{\perp}(\omega)&:= \int_{S^2} \rho(x,\omega,\omega_o)L_i(x,\omega)\,d\delta_{\omega_r}(\omega) \\
= \Lambda_{\omega_r}(\rho(x,\cdot,\omega_o)L_i(x,\cdot)) &= \rho(x,\omega_r,\omega_o) L_i(x,\omega_r).
\end{align}
The linear operator $\Lambda_{\omega_r}$ is not the BSDF itself, it is the operator $T$ when restricted to this specific specular vertex $x$, thus the BSDF cannot be a distribution even though Veach abuses notation (as is often done elsewhere) and calls it a distribution. What's really going on at a specular vertex $p$ is that the definition of $K$ at point $p$ is modified to $$(Kg)(p,\omega_o) = \int_{S^2}\rho(p,\omega,\omega_o)g(p,\omega)\,d\delta_{\omega_r}(\omega) = \rho(p,\omega_r,\omega_o)g(p,\omega_r).$$
Formally you could pretend that:
$$f(x,\omega,\omega_r) = \rho(x,\omega,\omega_o)\frac{d\delta_{\omega_r}(\omega)}{d\sigma_x^{\perp}(\omega)} =\rho(x,\omega,\omega_o)\delta_{\sigma_x^{\perp}}(\omega-\omega_r).$$
The only issue is that the function/Radon-Nikodym derivative doesn't exist because the Dirac measure $\delta_{\omega_r}$ is not absolutely continuous, but the notation is convenient enough that everyone abuses it.
Extension to Dirac Delta Light Sources
As we saw above the whole "distribution" business is tantamount to modifying your scattering operator at specific points. There are many ways to achieve this - you could tie the measure to your bsdf, then a specular vertex you would have a "differential bsdf-measure" $fd\delta$ while at a normal vertex it would be $fd\sigma_{x}^{\perp}$. The other way is of course to abuse notation and pretend that $\delta$ is the Radon-Nikodym derivative. The most theoretically sound, but least intuitive way, is to probably not touch the bsdf and just modify the scattering operator. Until this point we were lucky that we had an operator which to convert to our functional. It may not be readily obvious what should be done if we want a point or directional light source however.
As an example consider just the $0$-th term $L_e$ of the Neumann expansion with only a directional light source. To get the contribution from zero bounces away we integrate $L_e$ directly on the film. For a specific pixel on the film we may assume that the contribution is computed as $$\int_{A}\left(\int_{S^2}W_e(x,\omega)L_e(r(x,\omega),-\omega)\,d\sigma_{x}^{\perp}(\omega)\right)\,dx.$$
Here $A$ is the film area and $W_e$ is the sensitivity function. The aperture is implicitly taken into account by the ray-tracing function. But for simplicity assume there is no aperture, there are now two options: you either treat $L_e$ as an actual Dirac delta $L_e(y,\omega) = e(y,\omega)\delta(\omega-\omega_d)$ in which case you modify the integral to:
$$\int_{A}\left(\int_{S^2}W_e(x,\omega)e(r(x,\omega),-\omega)|n(x)\cdot\omega|\,d\delta_{-\omega_d}(\omega)\right)\,dx =\\= \int_{A}W_e(x,-\omega_d)e(r(x,-\omega_d),\omega_d)|n(x)\cdot\omega_d|\,dx.$$
You could however also require that $W_e$ is a Dirac delta, e.g. both w.r.t. position and direction. Once you do this the result will be a Dirac delta if you work with the abuse of notation approach (you will have one extra left over). This is not necessarily entirely wrong, if you assume that you really want a Dirac delta to mean that it carries infinite energy. As an example, if your sensor has both a positional (on the film) and directional Dirac delta, it would return the radiance of a normal light source when an importance ray hits it. However if that light source has a directional Dirac delta aligning with the one from the sensor, it is not that far fetched to say that you get an infinite contribution - and such an infinity that when integrated against a kernel would spit out the value at the kernel's specific location.
However, as far as I am aware we do not do this in computer graphics, e.g. when we have a point light source, directional light source, or even a laser, shining into the camera/film is not reproduced in the ray-tracers I am aware of. This means that they opted for treating this in a different way than the abuse of notation would suggest. E.g. if they have a specular light and then diffuse they treat it like a Dirac delta, however if they have a specular light and then a specular surface they do not use two Dirac deltas, instead they produce at most one Dirac delta, and that is only if the arguments of the two supposed Dirac deltas agree, so in this case $\int g(\omega)\delta(\omega-\omega_1)\delta(\omega-\omega_2)\,d\omega = g(\omega_1)\delta_{\omega_1=\omega_2} \ne g(\omega_1)\delta(\omega_1-\omega_2)$, note that in the middle expression the delta is a Kronecker delta and not a Dirac delta. The middle expression is consistent with how you do not visualize point/directional light sources on the film.
Conclusion
Ultimately how you decide to treat it is up to you, probably the most consistent way is using the abuse of notation, because then you can actually have multibounce specular contributions from a Dirac light source (e.g. "directional light" seen in a mirror), and it's also consistent as a limiting case of decreasing the solid angle emission of some light source and increasing its power - in the limit you get the directional Dirac delta, but the limiting emissions do not magically stop working after one bounce, so it's probably consistent to keep that for a Dirac delta light also. I think the main reason the Dirac delta nature of the light is not visualized on the film is because it is typically used as a cheap approximation of an actual physical light source, and the fact that it has infinite energy at a point would probably be perceived as an artefact on the film.
$$$$
Appendix: Some Details on the Dirac Delta
Here are some additional clarifications regarding the object (or I should rather say notational tool) that we call the Dirac delta. Essentially we wish for a function $\delta : \mathbb{R} \to \mathbb{R}\cup \{\pm\infty\}$ such that for any function $g$ we have $\int_{\mathbb{R}} g(t)\delta(t-x)\,dt = g(x)$. Unfortunately, such a function doesn't exist, however, there are objects (that are not functions) that behave in this way. I will list those roughly in order of most general to least general. The most general ones have the least assumptions on the class of functions $g$ for which they hold, while the least general ones have the most constraints on $g$.
Evaluation Functional $\delta_x$ (Constraint: Functions)
Let $V$ be the set of functions from some set $X$ to some set $Y$, i.e. $g:X\to Y$. Define $\delta_x(g) = g(x)$ - in this specific case we didn't assume anything about $g$ except that it is a function. If $Y$ is a vector space over a field $\mathbb{F}$, then $V$ is also a vector space over $\mathbb{F}$ with the function-function addition $(g+h)(x) := g(x)+h(x)$ and scalar-function multiplication $(\alpha\cdot g)(x) := \alpha\cdot g(x)$. Specifically if $Y=\mathbb{F}$ then $\delta_x : V \to \mathbb{F}$ is from the algebraic dual $V'$ of $V$, that is, the space of linear functions $f^*:V\to \mathbb{F}$ (also called functionals). Then $\delta_x (g) = g(x)$ is known as the evaluation functional (here we just assumed that $g$ are scalar functions). However you should notice that in this most general case we do not even have a notion of integration, i.e. we just have a notion of pointwise evaluation $\delta_x(g) = g(x)$ - the integral from our desired identity has been replaced by some map $\delta_x$.
Note also that $V$ being a set of functions is crucial. If I had taken $V$ to be $L_2(X)$ then I would be working with equivalence classes of functions and even the evaluation functional and the meaning of $g(x)$ would not be well-defined. One could of course try and pick out a representative from each equivalence class to which to apply the evaluation functional, but there is generally no well-defined way to do so for $L_2$. There is, however, a way to generalize this to some other spaces of equivalence classes, e.g. you could do so for Sobolev spaces $H^{q}(\mathcal{X})$ where $q\in(0,\infty)$ and $2q> \dim X$ (the dimension here is in the sense of manifold dimension).
Measure $\delta_x$ (Constraint: Measurable Functions)
If you have two measurable spaces $(X,\mathcal{X})$ and $(Y,\mathcal{Y})$ (those can be equal) and a function $g:X\to Y$ such that the preimage $g^{-1}$ of any set $B$ from $\mathcal{Y}$ is from $\mathcal{X}$: $B\in \mathcal{Y} \implies g^{-1}(B)\in \mathcal{X}$, then $g$ is termed a measurable function. If you additionally have a measure $\mu$ (e.g. the Lebesgue measure) then you can further limit $\mathcal{X}$ to be sets of Lebesgue-measurable sets, and we speak of $\mu$-measurable functions. For any $A\in \mathcal{X}$ you could now define the Dirac measure:
$$\delta_x(A) = 1_A(x) = \begin{cases} 1, & x\in A, \\ 0, & x\not\in A.\end{cases}$$
Then you can define the integral of $g$ w.r.t. $\delta_x$ as:
$$\int_{A} g(t)\,\delta_x(dt) = g(x)\delta_x(A)=
\begin{cases}
g(x), & x\in A, \\
0, &x\not\in A,
\end{cases}$$
Note that $0$ should be well-defined in $Y$ for this to make sense (all the practical examples you care about satisfy this however).
Now take $X$ to be some subset of $\mathbb{R}^n$ and $Y\subseteq \mathbb{R}$ then let $\mathcal{X}$ be the Lebesgue measurable subsets of $X$, and let $\mu:\mathcal{X}\to\mathbb{R}_+$ be the Lebesgue measure, then you can define integrals of the form $\int_A g(t) \,\mu(dt)$ for $A\in\mathcal{X}$. Each $\mu$-measurable function $f$ can then be interpreted as a density and it induces a measure $\mu_f(A) = \int_A f(t)\,\mu(dt)$. The Radon-Nikodym theorem goes in the opposite direction - it says that if you have two measures $\mu$ and $\nu$, then under some constraints (the measures being sigma-finite and $\nu$ being absolutely continuous w.r.t. $\mu$) there exists a $\mu$-measurable density function $f$ such that $\nu(A) = \int_A f(t)\,\mu(dt)$. Then $f$ is termed the Radon-Nikodym derivative of $\nu$ w.r.t. $\mu$: $f(t) = \frac{\nu(dt)}{\mu(dt)}$ and you can compute integrals w.r.t. $\nu$ by rewriting them as integrals w.r.t. $\mu$ (that is e.g. how you define $n$-volume integration on $n$-manifolds): $$\int_A g(t)\nu(dt) = \int_A g(t)f(t) \,\mu(dt).$$
In any case, using the above we see that if we have many compatible measures we can pick one common measure $\mu$ (e.g. the Lebesgue measure) and define all integrals w.r.t. it. It's what you do in practice, except when you encounter a Dirac delta, but then you modify your notation to look as if you're still using integration w.r.t. the Lebesgue measure for the Dirac delta. In our specific case it would be convenient if we could write
$$\int_{A} g(t)\,\delta_x(dt) = \int_A g(t)\delta_{\mu}(t-x)\,d\mu(t) = g(x)\delta_{x}(A).$$
And that would have been possible if it weren't for the fact that $\delta_x$ is not absolutely continuous w.r.t. $\mu$. Notably $\mu(\{x\}) = 0$ but $\delta_x(\{x\})=1$ which violates absolute continuity. This means that such a measurable density function $\delta_{\mu}(t-x)$ does not exist/the Radon-Nikodym derivative $\frac{d\delta_x}{d\mu}$ does not exist. But since it is convenient notation, one goes ahead and defines the formal symbolic expression by analogy to the cases where the Radon-Nikodym derivative exists:
$$\int_{A}f(t)\delta_{\mu}(t-x)\,\mu(dt) = \int_{A}f(t)\frac{\delta_x(dt)}{\mu(dt)}\,\mu(dt) := \int_A f(t)\,\delta_x(dt) = f(x)\delta_x(A).$$
The first two expressions are purely symbolic and equal to what's on the right by definition. Here $\mu$ is the Lebesgue measure or any other measure derived from the Lebesgue measure such as the solid angle or area measure used in the rendering equation. Note that you could try and define:
$$\delta_{\mu}(x-t) := \frac{\delta_x(dt)}{\mu(dt)}, \text{ but } \frac{\delta_x(dt)}{\mu(dt)} \text{ doesn't exist.}$$
So the $\delta_{\mu}$ can only show up in an integral, and then it's just symbolic notation that translates to what we defined it as. Notably the following expression (potentially corresponding to a specular light and bsdf) is not defined:
$$\int_{A} g(t)\delta_{\mu}(t-x)\delta_{\mu}(t-y)\,d\mu(dt) = \int_{A} g(t)\delta_{\mu}(t-y)\,\delta_x(dt) = ?,$$
simply because we never defined what $\int_{A} \delta_{\mu}(t-y)\,\delta_x(dt)$ means - remember that $\delta_{\mu}(t-x)$ doesn't actually exist and we only introduced notation for it w.r.t. an integral w.r.t. $\mu$ $\mu$. We could try and define a second "integral" of the above (remember that this is purely symbolic notation, so we're just mimicking what would have happened if the Radon-Nikodym derivative existed):
\begin{align}
I &=\int_B \int_{A} g(t)\delta_{\mu}(t-x)\delta_{\mu}(t-y)\,d\mu(dt)\,d\mu(dx) \\
&= \int_B\int_{A} g(t)\delta_{\mu}(t-y)\,\delta_x(dt)\,d\mu(dx) \\
&= \int_B g(x)\delta_{\mu}(x-y)1_A(x)\,\mu(dx) \\
&= \int_B g(x)1_A(x)\,\delta_y(dx) \\
&= g(y)1_A(y)1_B(y) = g(y)1_{A\cap B}(y).\end{align}
So we could claim that the result is a differential measure element $g(x)1_A(x)\,\delta_y(dx)$ that expects to be integrated. Later on we can argue what this should mean if we treat $dx$ as a finite measure set and we shrink it down to zero - then naturally you would expect this result to blow up to infinity.
However, whenever you have a specular light + specular bsdf, another way to model it is as follows:
\begin{align}
I &=\int_{A} g(t)\delta_{\mu}(t-x)\delta_{\delta_x}(t-y)\,d\mu(dt) \\
&= \int_A g(t)\delta_{\delta_x}(t-y)\,d\delta_x(dt) \\
&= \int_A g(t)\delta_{\delta_x}(t-y)\frac{\delta_y(dt)}{\delta_x(dt)}\,d\delta_x(dt) \\
&= \int_A g(t) \delta_y(dt) = g(y)1_A(y).
\end{align}
Now this corresponds to what is done in graphics $\delta_{\delta_x}$, but I think that the first approach is more consistent with what would happen in reality. I believe that what we do in graphics is motivated more from practical considerations than trying to model things realistically. For example if you were to take a disk light source with directional Dirac deltas for each point of it, and if you have complex ideal mirrors after it until the film, which say is also such a Dirac delta directional disk, you would have to track an infinite amount of specular rays. You could of course try to Monte Carlo this, but if you say that these carry infinite energy then missing even just one can result in an infinite error at the film.
In either case, since the Radon-Nikodym derivative doesn't actually exist both approaches amount to just changing the integration measure associated with the operator $K$ at specific points of the scene - so in theory you are free to pick whichever you like.
Distribution $\delta_x$ (Constraint: $C^{k}$ Functions)
You may require $\delta_{x}$ to be from the continuous dual $V^*$ (you need to have a topology to define the continuous dual unlike for the algebraic dual). If $V=C^{\infty}_c(X)$ then the continuos dual $V^*$ is the space of distributions. A regular distribution is one arising as an integral of a locally $L_1$ integrable function, e.g. $f^*(\phi) = \int f(t)\phi(t)\,dt$. You can multiply a distribution $D$ by another distribution $f^*$, if the latter is a regular distribution arising from $f\in C^{\infty}(X)$. So if your radiance and BSDF were infinitely many times differentiable you could make distributions out of those and multiply those distributions, similar to how you would deal with functions. The infinitely many times differentiability condition can be relaxed if you consider distributions of finite order. What cannot be relaxed however is the fact that the Dirac delta distribution is not a regular distribution - meaning that there exists no functions $\delta$ such that $\delta_x(\phi) = \int \delta(t-x) \phi(t)\,dt$. So once again such notation would be purely symbolic. Furthermore while distribution multiplication is defined between a regular distribution (arising from a function in $C^k$) and a non-regular one, the multiplication between two arbitrary distributions is not well-defined (e.g. between two Dirac delta distributions). As before you can try to define what it should be in your case depending on how you want to model your problem.
Reproducing Kernel $\delta_x$ (Constraint: Reproducing Kernel Hilbert Spaces)
I should note that the Dirac delta exists as a function in the sense $\int \delta(t-x) f(t)\,dt = f(x)$ if you consider a reproducing kernel Hilbert space. Then the Dirac delta function could be identified with the reproducing kernel for the space. For example bandlimited continuous functions have the sinc kernel as a reproducing kernel. Unfortunately you cannot guarantee that your functions are bandlimited and continuous in graphics - any hard edge can produce a non-bandlimited or discontinuous function in your radiance, or e.g. an edge in the BSDF can do the same (e.g. checkerboard texture).
The Dirac Delta as a Limit
You can study the following limit:
$$\lim_{\sigma\to 0}\int_{-\infty}^{\infty}f(t)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt = f(x).$$
As you can see it acts as a Dirac delta (for reasonable $f$). One may be tempted to move the limit inside the integral and state:
$$\delta(t-x) = \lim_{\sigma\to 0}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right),$$
but you actually cannot do so as the two are not equivalent:
$$\int_{-\infty}^{\infty}f(t)\lim_{\sigma\to 0}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt = 0
\ne f(x) = \lim_{\sigma\to 0}\int_{-\infty}^{\infty}f(t)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt.$$
However as before you could formally define the following symbolic expression:
$$\int_{-\infty}^{\infty}f(t)\delta(t-x)\,dt := \lim_{\sigma\to 0}\int_{-\infty}^{\infty}f(t)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt.$$
Then you could choose to define the product of two Dirac deltas as:
$$\int_{-\infty}^{\infty}\delta(t-x)\delta(t-y)\,dt := \lim_{\sigma_1,\sigma_2\to 0}\int_{-\infty}^{\infty}\frac{1}{2\pi|\sigma_1\sigma_2|}\exp\left(-\frac{(t-x)^2}{2\sigma_1^2}+\frac{(t-y)^2}{2\sigma^2_2}\right)\,dt.$$
My guess is that this integral should give you the following (I haven't checked it):
$$\lim_{\sigma\to 0}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right) = \begin{cases} \infty, & x=y, \\ 0, &x\ne y\end{cases}.$$
Based on the above and based on the idea of representing a Dirac delta light source as a limit of physical light sources where you increase the power while decreasing the solid angle or/and area, you could argue that e.g. a directional light source should show up as infinity if the film also has a directional Dirac delta importance/sensor function. So probably the first approach discussed in the section on $\delta_x$ as a measure is the physically more relevant one, since it would agree with this limiting process.
ifs), the mathematics are clearly defined, right? Looking forward for somebody who can clarify that. Besides, in nature there are no real dirac deltas / infinity values since there are no perfect mirrors; but I guess that is another topic. $\endgroup$