How non-abelian gauge coupling runs below confinement or QCD scale? I know the famous beta function of asymptotic free, but that seems describe the running coupling beyond confinement/QCD scale so that a perturbative analysis can apply. But how coupling runs below that scale? Any comment or references are greatly appreciated.
 A: Below that regime, we have the strongly coupled regime where perturbative approaches fail, due to the large value of coupling constant $\alpha_S$. The same is related to the QCD $\beta$ function via this relation. The behavior as a function of the energy scale looks roughly like this. Any perturbation expansion in this regime would give a divergent series, higher order terms contributing more than the lowest order. While that causes a lot of inconvenience (since we can't use our favorite and commonly employed tool, Feynman calculus), all is not lost - non-perturbative approaches are still applicable in this regime. e.g. Variational methods, based on minimization of the energy density, thermodynamic potentials etc. to understand the properties of the vacuum, are known to have been employed in toy models like the Gross-Neveu model. 
Another frequently used tool is to use effective hadronic models tailored for this context. Owing to quark confinement, if we have a bunch of quarks in this low-energy, strong coupling regime of QCD, they will in all probability be found as quark bound states instead of free quarks - i.e. mesons (${\bar q}q$) and baryons ($qqq$). Thus, taking them to be the effective degrees of freedom in this regime (called the hadronic regime, for obvious reasons), one can construct an effective theory of interacting baryons and mesons, which would approximate QCD at low-energies. The systematic recipe for these sort of things can be found in the Chiral Perturbation Theory ($\chi PT$) approach. (A more systematic introduction to this method can be found here.) While one can use this approach to get answers at various orders - leading order (LO), NLO, NNLO, sometimes even beyond, a sister approach is to get answers at LO, or NLO, and thereafter resort to phenomenology. i.e. you fit the free parameters of the LO chiral Lagrangian using some physical constraints, and thereafter use this Lagrangian to calculate observables of interest to you. 
A very popular and successful method often employed in the past two decades, is the Walecka model, wherein you write down a similar effective Lagrangian for this many-body problem. (An old review of the method can be found here.) The original idea was to introduce a less massive scalar meson (like $\sigma$ meson) to provide a long-range attractive interaction, and a more massive vector meson (like $\omega$ meson) to provide a short-range repulsive interaction. The same strategy was subsequently generalized to include chiral invariance as well, akin to the $\chi PT$ approach, in the so-called chiral $\sigma-\omega$ models.
So, long story cut short - mankind has evolved specific methods for tackling this strong coupling regime and the rich physics residing here.
