with \(L = \ln (Q^2/\Lambda _{\rm QCD}^2)\) . Here the integration constant for solving (REF ) is chosen by
\((\beta _1/\beta _0^2) \ln (\beta _0)\)  [1]}. The expansion coefficients of the
\(\beta \) -function to N\(^3\) LO were calculated in [2]}.
The flavor matching conditions were given in [3]}.
The expansion coefficients of the renormalized mass were given in [4]}, [5]}.
The constant and \(O(\varepsilon )\)  parts of the massive unrenormalized operator matrix elements at \(O(a_s^k)\) 
are denoted by \(a_{ij}^{(k)}\)  and \(\bar{a}_{ij}^{(k)}\) , respectively, cf. [6]}, [7]}, [8]}, [9]}.
