normalized Laplacian of Gaussian Laplacian of Gaussian formula for 2d case is
$$\operatorname{LoG}(x,y) = \frac{1}{\pi\sigma^4}\left(\frac{x^2+y^2}{2\sigma^2} - 1\right)e^{-\frac{x^2+y^2}{2\sigma^2}},$$ in scale-space related processing of digital images, to make the Laplacian of Gaussian operator invariant to scales, it is always said to normalize $LoG$ by multiplying $\sigma^2$, that is
$$\operatorname{LoG}_\text{normalized}(x,y) = \sigma^2\cdot \operatorname{LoG}(x,y) = \frac{1}{\pi\sigma^2}\left(\frac{x^2+y^2}{2\sigma^2} - 1\right)e^{-\frac{x^2+y^2}{2\sigma^2}}.$$ I wonder why multiply by $\sigma^2$ not $\sigma^4$ or anything else?
UPDATE
Thanks to comments from @achille. From the perspective of dimensional analysis, in the Laplacian of Gaussian operator
$$LoG(x,y,\sigma)=\frac{\partial^2g}{\partial x^2} +\frac{\partial^2g}{\partial y^2}$$, I think $x,y$ are variables with dimension $L$, $\sigma$ is a parameter with dimension $L$. But what about $g$? Since $g$ is a function of $x,y,\sigma$,
$$g(x,y,\sigma)=\frac{1}{2\pi \sigma^2}exp(-\frac{x^2+y^2}{2\sigma^2})$$,
and $x,y,\sigma$ are of the same dimension $L$, so I guess in $g$, the term $exp(-\frac{x^2+y^2}{2\sigma^2})$ is dimensionless, isn't it? And the term $\frac{1}{2\pi \sigma^2}$ is of dimension $L^{-2}$, right? So $g$ is actually of dimension $L^{-2}$, isn't it?
Now come back to $LoG$, it should have dimension $L^{-4}$?
UPDATE 2
Laplacian operator is 
$$\nabla^2 = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$$,
and a Gaussian function's scale is $\sigma$, right? If I apply $\nabla^2$ on a Gaussian function $g(x,y,\sigma)$, what is difference of applying the dimensionless $\sigma^2*\nabla^2$?
 A: First, let me try to give you some intuition of why you have to normalize by scale at all. As you go from finer to coarser scales you blur the image. That makes the intensity surface more and more smooth. That, in turn, means that the amplitude of image derivatives gets smaller as you go up the scale volume.  This is a problem for finding interest points, because you are looking for local extrema over scale. Without normalization you will always get the maximum at the finest scale and the minimum at the coarsest scale, and that's not what you want.
So, image derivatives are attenuated as $\sigma$ increases. In fact, the derivatives decrease exponentially as a function of $\sigma$. To compensate for that you have to normalize them by multiplying the $n$-th derivative by $\sigma^n$. Since the LoG is a combination of second derivatives, you have to multiply it by $\sigma^2$. 
You can find the derivation and a better explanation of this in this paper by Toni Lindeberg.
A: EDIT As far as I can (now) tell, your problem is unrelated to dimensional analysis, so yesterday's input (below) is obsolete. Now, your LoG function is the Laplacian of
$$
G(x,y;\sigma) = \frac{1}{2\pi\sigma^2} {\rm e}^{-\frac{x^2+y^2}{2\sigma^2}} , 
$$
whence
$$
{\rm LoG}(x,y;\sigma)
=
\Delta_{(x,y)}G(x,y;\sigma)
=
\frac{\partial^2 G(x,y;\sigma)}{\partial x^2} + \frac{\partial^2 G(x,y;\sigma)}{\partial y^2}
=
\frac{1}{\pi\sigma^4} \left(\frac{x^2+y^2}{2\sigma^2}-1\right){\rm e}^{-\frac{x^2+y^2}{2\sigma^2}} .
$$
Multiplying the latter by $\sigma^2$ is equivalent to multiplying the former by the same number, so in fact you're working with
$$
G_{\rm norm}(x,y;\sigma) = \frac{1}{2\pi} {\rm e}^{-\frac{x^2+y^2}{2\sigma^2}}
$$
and
$$
{\rm LoG}_{\rm norm}(x,y;\sigma)
=
\Delta_{(x,y)}G_{\rm norm}(x,y;\sigma)
=
\frac{1}{\pi\sigma^2} \left(\frac{x^2+y^2}{2\sigma^2}-1\right){\rm e}^{-\frac{x^2+y^2}{2\sigma^2}} .
$$
Note that both $G$ and $G_{\rm norm}$ are rescalings of the standard (i.e., unit $L^1-$norm) Gaussian. Knowing nothing about image processing, the one interesting thing that strikes me about ${\rm LoG}_{\rm norm}$ is the following identity:
$$
\Delta_{(x,y)}G_{\rm norm}(\sigma x,\sigma y;\sigma) = \Delta_{(x,y)}G_{\rm norm}(x,y;1) .
$$
The non-normalized version unsurprisingly has a factor of $\sigma^2$ in front of the LHS. Make what you want of this i.t.o. image processing, I'm afraid I can't help without a clear mathematical objective.
OBSOLETE There's a fundamental misunderstanding here related to what you write: 'I guess in g, the term ${\rm exp}\left(−\frac{x^2+y^2}{2σ^2}\right)$ is dimensionless [...]' No, this does not need be the case: the unit in front of the exponential - the one ($1$) multiplying it - is the one carrying the units of $g$.
To close with a simple example, recall that a body traveling with 'constant velocity equal to one' travels a distance $s(t)=t$ after time $t$. Does that mean that distance $s$ has the units of time $t$? Of course not; it means that $s(t)=vt$, with the velocity $v$ here being equal to one in the chosen system of units. Change units and $v$ will become something else: for example, if $v=1{\rm m/sec}$ and $t$ is measured in seconds, then $s(t)=t$ as long as $s$ is measured in meters: after $t$ seconds, the body has traveled $t$ meters. Change from meters to millimeters, so that $v=1000{\rm mm/min}$ and $s(t)=1000t$: after $t$ seconds, the body has traveled $1000t$ millimeters. To students learning the concept of dimensionality, it tends to be clearer that $1000$ carries the units (${\rm mm/sec}$) in the latter problem than that the unit does the same in the former one - that's because they forget that the unit is even there.
