Abstract. Dyson-Schwinger equations have been successful in determining the correlation functions in Yang-Mills theory in the Landau gauge, in the infrared regime. We argue that similar results can be obtained, in a technically simpler way, with Callan-Symanzik renormalization group equations. We present generalizations of the infrared safe renormalization scheme proposed by Tissier and Wschebor in 2011, and show how the renormalization scheme dependence can be used to improve the matching to the existing lattice data for the gluon and ghost propagators.Historically, the first semi-analytical method that has made it possible to access the deep infrared (IR) regime of Landau gauge Yang-Mills theory, have been Dyson-Schwinger equations. The first solutions of these equations that have been found show a scaling behavior of the propagators in the IR [1]. Several years later, another type of solutions of the same equations was discovered [2], with a massive gluon propagator and a finitely enhanced ghost propagator in the IR (decoupling solutions). To date, Dyson-Schwinger equations are still the main tool for the (semi-)analytical exploration of the IR regime of Yang-Mills theory.In a parallel development that actually initiated much earlier with Gribov's observation of the existence of gauge copies in his famous 1978 paper [3] and was later worked out in great detail by Zwanziger [4], the theoretical foundations of IR Yang-Mills theory in the Landau gauge were laid. While this so-called Gribov-Zwanziger scenario seemed to confirm the existence of the scaling solutions, a more recent refinement [5] favors the decoupling type of solutions. Finally, the results of simulations of Yang-Mills theory in the Landau gauge on huge lattices [6] have been interpreted by most workers in the field as confirming the realization of the decoupling type of solutions in three and four space-time dimensions.In this contribution, we will present a different technique for the exploration of the IR regime of Yang-Mills theory. We intend to reproduce the results of lattice simulations in the Landau gauge that restrict the gauge field configurations to the (first) Gribov region, but make no attempt to reach the fundamental modular region. The restriction to the Gribov region implies the breaking of BRST invariance in the continuum formulation [4]. The most important consequence of the broken BRST invariance for the IR regime is the appearance of a mass term for the gluon field [5]. Indeed, the addition of a gluonic mass term to the Yang-Mills action has been shown to reproduce all the solutions (two types of scaling solutions and the decoupling solution) found before with the help of DysonSchwinger equations when solving the Callan-Symanzik equations for this theory in the IR regime in an epsilon expansion around the upper critical space-time dimension [7]. This analysis also shows that a