What is the fundamental reasoning for anti aliasing using multiple random samples within a pixel?
From a signal processing point of view, you're sampling a continuous-domain signal, and you need to filter it to get rid of frequencies beyond the Nyquist limit. It's that filtering that leads to integrating over the pixel area—or more generally, integrating over the support of your antialiasing kernel (which need not be a box).

Consider your rendering function that takes a sample point $x, y$ in screen space and returns back the color found at that point. (Let's ignore any issues of random sampling for the moment, and assume it returns a "perfect" rendered color for that specific point.) This function effectively defines a 2D continuous-domain signal. Or to put it another way, it defines an infinite-resolution image, since nothing prevents this function from having features at arbitrarily small scales. In terms of frequency domain: the function is not band-limited; it can include components of arbitrarily high spatial frequencies.

Now you want to convert it down to a finite number of pixels. Much like digitizing an audio signal, when you sample it, you'll get aliasing unless you first eliminate frequencies beyond the Nyquist limit imposed by the sampling rate. In other words, you have to get rid of features smaller than the pixel grid. To do this, you apply a low-pass filter. The ideal low-pass filter is the sinc function, but for various reasons of practicality we use other filters (that don't perfectly eliminate frequencies beyond the Nyquist limit, but they at least attenuate them).

Low-pass filtering is done by convolution. If $f(x, y)$ is the rendering function and $k(x, y)$ is a filter kernel, then we can define the low-pass filtered image as:

$$f_\text{filtered}(x, y) = \iint f(x', y') \, k(x' - x, y'- y) \, dx' \, dy'$$

Then the image is safe to sample, so the final pixel values are gotten by just evaluating $f_\text{filtered}$ at the pixel coordinates.

If $k$ is a box filter, which looks like $k = 1$ within the pixel box and $k = 0$ elsewhere, then this simplifies to just integrating $f$ over the pixel box. But as noted, the box filter isn't so great and there are better choices such as tent, bicubic, and Gaussian filters.

Anyway, now that we have an integral, we can use Monte Carlo for it and combine it with any other integrals we might like to do—for lighting, motion blur, and so on. We can even apply importance sampling to the $k$ factor in the integral, by generating samples for each pixel that are distributed around the pixel center according to $k$.