  We derive the multi-fractal scaling of probability distributions of
multi-particle configurations for the binary reaction-diffusion system $A+A \to
\emptyset$ in $d \leq 2$ and for the ternary system $3A \to \emptyset$ in
$d=1$. For the binary reaction we find that the probability $P_{t}(N, \Delta
V)$ of finding $N$ particles in a fixed volume element $\Delta V$ at time $t$
decays in the limit of large time as $(\frac{\ln t}{t})^{N}(\ln
t)^{-\frac{N(N-1)}{2}}$ for $d=2$ and
$t^{-Nd/2}t^{-\frac{N(N-1)\epsilon}{4}+\mathcal{O}(\ep^2)}$ for $d<2$. Here
$\ep=2-d$. For the ternary reaction in one dimension we find that
$P_{t}(N,\Delta V) \sim (\frac{\ln t}{t})^{N/2}(\ln
t)^{-\frac{N(N-1)(N-2)}{6}}$. The principal tool of our study is the dynamical
renormalization group. We compare predictions of $\ep$-expansions for
$P_{t}(N,\Delta V)$ for binary reaction in one dimension against exact known
results. We conclude that the $\ep$-corrections of order two and higher are
absent in the above answer for $P_{t}(N, \Delta V)$ for $N=1,2,3,4$.
Furthermore we conjecture the absence of $\ep^2$-corrections for all values of
$N$.
