In this chapter we derive the analytical formalism allowing for
the all order \(\varepsilon \) -expansions for the integrals up to two denominators (i.e. 
the basic angular integrals) in terms of the single and double
Nielsen polylogarithms [1]}, [2]}.
As it was shown
in Sec. REF  for arbitrary dimension \(D\) 
the double massive angular integral can be transformed to a finite sum of
single-massive integrals based on the reduction formalism proposed by us. Therefore,
\(\varepsilon \) -expansion of the double massive integral for any \(j,l\in \mathbb {Z}\)  can be generated
by use of the \(\varepsilon \) -expansion of the single massive integral.
Using two mass splitting (see Sec. REF )
and recursion relations (see Sec. REF )
any angular integral \(I_{j,l}(v_{12},v_{11},v_{22};\varepsilon )\) 
with \(j,l\in \mathbb {Z}\)  can be expressed in terms of the basic angular integrals
mentioned above.
Note, that these results have never been discussed or derived before in literature.
Only expansion up to a few first orders in \(\varepsilon \)  were known before
(see, e.g. results in Refs. [3]}, [4]}, [5]}).
Therefore, we think that our findings will be useful as input for the calculation
of QCD processes.
