What is exactly a quantum fluctuation? According to the quantum fluctuation by Heisenberg Uncertainty Principle, matters and antimatters are created and disappearing by colliding in an almost complete vacuumed area of space near the event horizon. However, isn't this cannot happen because the mass cannot be created nor destroyed by the law of conservation of mass?
 A: Quantum fluctuation is the phenomenon of observing random outcomes from an experiment due to the quantum uncertainty principle. More rigorously, an observable in a given instant of time $t$ is described in quantum mechanics by a Hermitian operator $\mathcal O(t)$ and the theory cannot predict the precise value for this quantity after a given measurement. It only predicts the probability of various outcomes by computing the trace:
$$
P_{\mathcal O\rightarrow o}=tr(\rho_{\psi}|o(t)\rangle\langle o(t)|)
$$
where $\rho_\psi = |\psi\rangle\langle\psi|$ is the state of the system, an object that is determined by the preparation of the system, and $|o(t)\rangle$ is a vector that satisfies:
$$
\mathcal O(t)|o(t)\rangle= \mathcal o(t)|o(t)\rangle
$$
where $o(t)$ is a possible value that $\mathcal O (t)$ can take at time $t$. Here I am using the Heisenberg prescription for the time evolution.
In some very special situations turn outs that the probability of a given outcome is one or close to one, such that there is no "fluctuation" taking place in the outcomes. Now, if the probability is far from one we expect that we will see fluctuations in our outcomes if we manege to reproduce the experiment $N$ times.
The source for the fluctuations are the non-commutative properties
of the hermitian operators that describes the observables. For example, the non-commutativity between the observable we are measuring and the state associated to the preparation of the experiment
$$
[\rho_{\psi}, \mathcal O]\neq 0
$$
which can be understood as some kind of incompatibility between the preparation of the experiment and the quantity we are measuring, like preparing a spin particle with a given value of $S_z=+1/2$ and measuring the $S_{x}$ component.
Now, this does not enter in conflict with energy conservation because the value of energy is only well defined in very special cases. In quantum mechanics this means that the hermitian operator $H(t)$ that is associated with energy, the Hamiltonian, may not commute with the state $\rho_{\psi}$. The conservation of energy in quantum mechanics takes a very precise form and is given by the following equation
$$
H(t)=H(0)
$$
for every instant of time $t$. This follows from the fact that the equation that describes the time evolution of the observables is given by:
$$
\frac{d\mathcal O}{dt}=i[H(t),\mathcal O(t)]
$$
so for the case where $\mathcal O(t)=H(t)$ the RHS vanishes. This is valid for any closed system.
In Quantum Field Theory, the Hamiltonian does not commute with certain local operator, or operators that are too narrow in space-time, like the value of the charge distribution or energy distribution at short times. This will imply that such quantities will fluctuate both because they does not commute with the RHS of the equation that determines the time evolution and because the vacuum state $\rho_0$, a state with zero
energy, does not commute with such quantities.
A: Another way to think about quantum fluctuations is to say that "virtual particles" are field amplitudes appearing in important integrals, with contributions that don't obey $E^2=p^2+m^2$ aka $\omega^2=\vec{k}^2+m^2$. I'll summarize the relevant points from Chapters I.3-I.4 of Quantum Field Theory in a Nutshell, focusing on the simplest special case, but the idea is qualitatively the same for other interactions.
Work throughout in $+---$ with $c=\hbar=1$. Consider a mass-$m$ real scalar boson field $\varphi$ with current $J$, of Lagrangian density $\mathcal{L}=\frac12(\partial_\mu\varphi\partial^\mu\varphi-m^2\varphi^2)+J\varphi$ or, adding a total derivative without changing the physics, $\mathcal{L}=-\frac12\varphi(\partial_\mu\partial^\mu+m^2)\varphi+J\varphi$. The path integral$$Z(J):=\int\mathcal{D}\varphi\exp(i\int d^4 x\mathcal{L})=Z(J=0)\exp(iW(J))$$where$$W(J):=-\frac12\iint d^4xd^4yJ(x)D(x-y)J(y),$$with the propagator $D(x-y)$ solving $(\partial_\mu\partial^\mu+m^2)D(z)=-\delta^{(4)}(z)$. We can easily write $D$ as a Fourier transform, so$$W(J):=-\frac12\iiint d^4xd^4y\frac{d^4k}{(2\pi)^4}\frac{J(x)\exp (ik(x-y))J(y)}{k_\mu k^\mu-m^2+i\epsilon}.$$There are both mathematical and physical motives for the $i\epsilon$ term in the denominator. The former is to take care with contour integration; the denominator can't vanish for $\epsilon\in\Bbb R\setminus\{0\}$, but without this term we'd have zero denominator on-shell (i.e. when $k_\mu k^\mu=m^2$). The latter is to allow an arbitrarily small amount of damping in the physics, which we take to $0$ when doing so extracts a finite calculation of something. As an easy warm-up to the field-theoretic case above, note the solutions of $(\partial_t^2+\gamma\partial_t+\omega_0^2)y=\delta(t)$ include $\frac{1}{2\pi}\int\frac{\exp(i\omega t)d\omega}{\omega_0^2-\omega^2+i\gamma\omega}.$
Anyway, we can get an interaction between point sources at $\vec{x}_1\,\vec{x}_2$ from the cross terms for $J(x)=\sum_{a=1}^2\delta^{(3)}(\vec{x}-\vec{x}_a)$, viz.$$W_\text{int}(J):=-\iiint d^4xd^4y\frac{d^4k}{(2\pi)^4}\frac{\delta^{(3)}(\vec{x}-\vec{x}_1)\exp (ik(x-y))\delta^{(3)}(\vec{y}-\vec{x}_2)}{k_\mu k^\mu-m^2+i\epsilon}\\=\iint dx^0\frac{d^3\vec{k}}{(2\pi)^3}\frac{\exp i\vec{k}\cdot(\vec{x}-\vec{y})}{k_i k^i+m^2},$$where we can finally drop the $\epsilon$ from the denominator. Over a finite time $T$, the identification $\exp iW=\exp -i\int dx^0 E$ lets us identify the potential energy,$$E=-\int \frac{d^3\vec{k}}{(2\pi)^3}\frac{\exp i\vec{k}\cdot(\vec{x}-\vec{y})}{k_i k^i+m^2}=-\frac{1}{4\pi r}\exp(-mr),\,r:=|\vec{x}_1-\vec{x}_2|.$$This is the Yukawa potential; it reduces to a familiar-from-EM/gravity $-\frac{1}{4\pi r}$ potential, with force $-\frac{1}{4\pi r^2}$, if $m=0$. (Of course, those aren't due to scalar bosons, but that's not important right now.) But the exact formula as a function of $r$ isn't the issue; what matters is that $E$ ($W$) is an integral over $\vec{k}$ ($k$), not just over the "shell" satisfying the dispersion relation $k_\mu k^\mu=m^2$. We don't even need to keep $k_\mu k^\mu$ positive, just as in quantum tunneling you don't need a positive kinetic energy. However, as with the case of tunneling, the further we move into the classically forbidden region the less contribution the scenario makes to what happens on average.
But as @Nogueira noted, an off-shell analysis is unnecessary. We can instead obtain momentum-space transition amplitudes, e.g. Eq. (7) here gives a matrix element $\langle\vec{p}_1^\prime\vec{p}_2^\prime|M_A+M_B|\vec{p}_1\vec{p}_2\rangle=\frac{g^2}{s-m^2}$ in terms of a Mandelstam variable. We can also not only discuss a field's quantum fluctuations on-shell, but can also distinguish them from thermal fluctuations. A Klein-Gordon field $\phi_t(\vec{x})$ of Fourier transform $\tilde{\varphi}_t(\vec{k})$ gives fluctuations of probability density $\exp-\int\frac{d^3\vec{k}}{(2\pi)^3}\tilde{\varphi}_t^\ast(\vec{k})\omega_\vec{k}\tilde{\varphi}_t(\vec{k})$. (See here for the thermal counterpart.) 
