Abstract: There is a common belief that the main uncertainties in the theoretical analysis of neutrinoless double beta (0 nu beta beta) decay originate from the nuclear matrix elements. Here, we uncover another previously overlooked source of potentially large uncertainties stemming from nonperturbative QCD effects. Recently perturbative QCD corrections have been calculated for all dimension 6 and 9 effective operators describing 0 nu beta beta-decay and their importance for a reliable treatment of 0 nu beta beta-decay has been demonstrated. However, these perturbative results are valid at energy scales above similar to 1 GeV, while the typical 0 nu beta beta scale is about similar to 100 MeV. In view of this fact we examine the possibility of extrapolating the perturbative results towards sub-GeV nonperturbative scales on the basis of the QCD coupling constant “freezing” behavior using background perturbation theory. Our analysis suggests that such an infrared extrapolation does modify the perturbative results for both short-range and long-range mechanisms of 0 nu beta beta-decay in general only moderately. We also discuss that the tensor circle times tensor effective operator cannot appear alone in the low energy limit of any renormalizable high-scale model and then demonstrate that all five linearly independent combinations of the scalar and tensor operators, which can appear in renormalizable models, are infrared stable.