Definition of symmetrically ordered operator for multi-mode case? As I know, Wigner function is useful for evaluating the expectation value of an operator. But first you have to write it in a symmetrically ordered form. For example:
$$a^\dagger a = \frac{a^\dagger a + a a^\dagger -1}{2}$$
For single mode case where there is  only one pair of creation and destroy operator the symmetrically ordered operator is defined. But for multi-mode case,how is it defined? For example, how would we write 
$$a_1^\dagger a_1 a_2^\dagger a_2$$
in a symmetrically ordered form (such that we could easily evaluate its expectation value using Wigner function)?
 A: OK, I assume you are comfortable with the rules of Weyl symmetrization where you may parlay the $[x,p]=i\hbar$ commutation relation to the $[a,a^\dagger]=1$ one... the combinatorics is identical provided you keep track of the is and the ħs, etc...
So your surmise is sound. Since different modes commute with each other,
they don't know about each other, and you Weyl-symmetrize each mode factor separately. 
So, e.g., for your example, 
$$
a_1^\dagger a_1    a_2^\dagger a_2=   a_2^\dagger a_2 a_1^\dagger a_1      = \frac{a_1^\dagger a_1 + a_1 a_1^\dagger -1}{2}  ~ \frac{a_2^\dagger a_2 + a_2 a_2^\dagger -1}{2},
$$
etc.
All you need to recall is that each mode is in a separate (Fock) space of a tensor product, so all operations factor out into a tensor product of Weyl-symmetrized factors.
This is the reason all multimode/higher-dim phase space generalizations of all these distributions functions are essentially trivial, and most subtleties are routinely illustrated through just one mode. 

Edit in response to comment on entangled states. One of the co-inventors of the industry, Groenewold, in his monumental 1946 paper, Section 5.06 on p 459, details exactly how to handle entangled states--in his case for the EPR state. The entanglement and symmetrization is transparent at the level of phase-space parameters (Weyl symbols): the quantum operators in the Wigner map are still oblivious of different modes. What connects/entangles them, indirectly, are the symmetrized δ-function kernels involved, even though this is a can of worms that even stressed Bell's thinking. The clearest "modern" paper on the subject is Johansen 1997, which, through its factorized Wigner function and changed +/- coordinates, reassures you you never have to bother with the quantum operators: the entangling is all in the Wigner function and phase-space, instead! (Illustration: 351884/66086.)
A: A quasiprobability distribution depends on the symbol/operator ordering prescription. E.g.: 


*

*For Weyl/symmetric ordering of Hermitian operators, one uses the Wigner quasiprobability distribution.

*For Wick/normal ordering of creation & annihilation operators, one uses the Glauber–Sudarshan $P$ representation.

*For anti-normal ordering of creation & annihilation operators, one uses the Husimi $Q$ representation.
A: Symmetrically order expansion of the ladder operator is written as follows;
a1b1=1/2(a1b1 +b1a1)= a1b1+1/2
a=creation operator b= anihilatinoperator, also
a1b1a2b2=1/2(a1b1 +b1a1)1/2(a2b2+b2a2)= (a1b1+1/2)(a2b2+1/2)=a1b1a2b2+1/2(a1b1)+1/2(a2b2)+1/4
A: Symmetrically order expansion of the ladder operator is written as follows;
$$a_1 b_1 =  (a_1 b_1 + b_1 a_1)/2= a_1 b_1 + 1/2$$
where $a$ is the creation operator and $b$ is the annihilation operator, also
$$a_1 b_1 a_2 b_2 = \frac{1}{2} (a_1 b_1 + b_1 a_1) \frac{1}{2}(a_2 b_2+b_2 a_2) = (a_1 b_1+ \frac{1}{2})(a_2 b_2+ \frac{1}{2}) = a_1 b_1 a_2b_2+  \frac{1}{2}(a_1 b_1)+\frac{1}{2}(a_2 b_2) + \frac{1}{4}$$
