Let us suppose I am taking some measurements whose outcomes are distributed according to a gaussian probability density distribution centered at $0$. For each measurement result falling in the range $\Delta x$ centered around result value $ x$ I then add 1 to the corresponding "placeholder". Of course, the resulting probability density plot (i.e. the plot of each placeholder count/[$\Delta x \ \times $ total number of measurements]) will be very well approximated (well, for reasonably small $\Delta x$) by a gaussian distribution, whose integral from $-\infty$ to $+\infty \, \,$ (i.e. total probability better be 1 ...) $$ \int_{-\infty} ^{\infty} \frac{e^{-\frac{x^2}{2\sigma^2}}}{\sqrt{2\pi} \,\,\sigma}dx\ = 1 $$ the above can be easily solved analytically, and the most naive numerical approximation of the above integral would actually correspond to directly summing together the probabilities corresponding to each placeholder. Defining $$ f(x) = \frac{counts \,\, at \,\, x}{total \,\, counts \,\, \times \Delta x} $$ said direct sum is $$ \sum_{-\infty}^{+\infty} \ f(x) \ \Delta x \ = \ 1 $$ and so far it has been quite trivial. But then I add a second identical measurement setup, and each time I take a measurement on the first setup I also take a measurement on the second one. Wishing to compute the total probability for the fraction of measurements corresponding to exactly the same placeholder for said pair of "simultaneous" measurements, I simply multiply together the two corresponding probabilities (each one of them being $f(x) \ \Delta \ x \ $) and add up: $$ \sum_{-\infty}^{+\infty} \ (f(x) \ \Delta x)^2 \ = \ ? $$ said $?$ value will certainly be $<1$, and it is easy to compute it numerically, adding up values in the above sum, while starting from negative $x$ values at which the corresponding placholder is practically empty, etc., etc. (when for the first sum I numerically compute $= 1$, for the second sum I obtain $? \ = \ 0.0282 ...$). Finally, I am getting to my question: naively, I would be tempted to approximate the above $?$ sum, whereby $f(x)$ is a gaussian probability density distribution, with $$ \int_{-\infty} ^{\infty} \Bigg( \frac{e^{-\frac{x^2}{2\sigma^2}}}{\sqrt{2\pi} \,\,\sigma}dx\ \Bigg)^{2} $$ but I am running into diffulties with that $(dx)^2$ now replacing the usual $dx$. Not being very knowledgeable about the fundamentals of the underlying measure theory, how can such an integral be approached ? but first of all, does it even make any sense at all ? and finally, is there any way to compute said $ \ 0.0282 ... \ $ value also analytically ? I have tried to play around by inserting the Dirac $\delta$ function, but got nowhere (the result shall be dimensionless, etc. etc.).
CLARIFICATION
I shall clarify that said $ \ 0.0282 ... \ $ value was obtained by numerically summing from $ x=-10$ to $x=+10$, and at discrete intervals of width $0.1$, for the particular case $\sigma = 1$. So that for example at $x=0.5$ the corresponding joint probability value would be $\ p_j(0.5) = (f(x) \ \Delta x)^2$ = 0.001239...
It is further interesting to remark how said $ \ 0.0282 ... \ $ value of the total probability for simultaneous identical measurements would then result multiplied by $k$ would $\sigma$ be divided by a value $k$. The searched for analytic expression would thus need to include an inverse proportionality wrt to $\sigma$.
I believe that I have now found an answer. By reasoning deeper about the actual physical meaning of the mentioned joint measurements, it would appear intuitive that by decreasing the size, $\Delta X$, of the "bin" at each placholder (and thus correspondingly increasing the total number of placeholders) the probability that a joint measurement results in two values falling in the same bin shall also decrease. Indeed, after running few more numerical computations, it finally appears that the correct value for said "same values" total probability is:
$$ 0.2820947917738780 \ \times \ \frac {\Delta X}{\sigma} $$ So, coming back to the original question, how could we approximate that with an integral expression ?
We may observe that the above sum should be reasonably approximated by the volume defined as the product of $\Delta X$ and the area under the function $\Big( \frac{1}{\sqrt{2\pi} \,\,\sigma} \ e^{-\frac{x^2}{2\sigma^2}}\Big) ^2 $, thus: $$ \sum_{-\infty}^{+\infty} \ (f(x) \ \Delta X)^2 \approx \int_{-\infty} ^{\infty} \Bigg( \frac{e^{-\frac{x^2}{2\sigma^2}}}{\sqrt{2\pi} \,\,\sigma}\Bigg)^2 \Delta X \ dx = \frac{1}{2\sqrt{\pi}} \frac{\Delta X}{\sigma} = 0.2820947917738780 \ \frac {\Delta X}{\sigma} $$ quite a reasonable approximation indeed (after all, for the numerical computations I was sampling out of gaussian distributions).