Renormalisation group evolution for the $\Delta S = 1$ effective Hamiltonian with $N_f=2+1$

We discuss the renormalisation group (RG) evolution for the $\Delta S = 1$ operators in unquenched QCD with $N_f = 3$ ($m_u=m_d=m_s$) or, more generally, $N_f = 2+1$ ($m_u=m_d \ne m_s$) flavors. In particular, we focus on the specific problem of how to treat the singularities which show up only for $N_f=3$ or $N_f = 2+1$ in the original solution of Buras {\it et al.} for the RG evolution matrix at next-to-leading order. On top of Buras {\it et al.}'s original treatment, we use a new method of analytic continuation to obtain the correct solution in this case. It is free of singularities and can therefore be used in numerical analysis of data sets calculated in lattice QCD.