How is $\Lambda_{\textrm{QCD}}$ relevant in the non-perturbative regime? The famous $\Lambda_{\textrm{QCD}}$ parameter enters through the one-loop running of the QCD coupling, through a relation similar to the following:
$$\alpha_S(Q^2)=\frac{\alpha_S(Q^2_0)}{1+b\ln(Q^2/Q^2_0)}\equiv\frac{\alpha_S(Q^2_0)}{b\ln(Q^2/\Lambda^2_\textrm{QCD})}$$
My question is simple: how does this equation, and thus any definition of $\Lambda_{\textrm{QCD}}$, have anything to do with QCD in the non-perturbative regime, where $\alpha_S>1$ and thus these equations break down? Here $Q$ is strictly just an arbitrary renormalization scale, but it could also be an energy scale in a particular process we are considering.
I know the conformal/trace anomaly in QCD is given be $T^\mu_\mu\sim \beta(\alpha_S)F^2$, where $\beta(\alpha_S)$ is the beta function. But is this an all-orders result? (i.e. are all the orders of $\alpha_S$ resummed correctly in $\beta(\alpha_S)$?) 
Also, I know that the mass of a hadronic state is given by the trace anomaly $\langle P|T|P\rangle\sim M^2$, but surely this cannot mean that $M^2\sim\beta(\alpha_S)$ because the beta-function is scheme dependent, whereas the mass of a hadron is totally physical.
But none of the two previous two paragraphs say anything about what happens when $Q\sim\Lambda_{\textrm{QCD}}$ and $\alpha_S\sim 1$. If this truly is how the $\Lambda_{\textrm{QCD}}$ parameter enters into life, I don't see how it isn't just a mere artifact of the logarithms which appear in the perturbative regime. It would then seem possible for its relevance to disappear if we somehow calculated the beta-function to all-loops.
 A: You are quite correct when you say that $\Lambda_{QCD}$ may be an artifact of perturbation theory. This actually is the current interpretation of $\Lambda_{QCD}$, based on a few observations.
$$
$$
Why we believe that $\boldsymbol{\Lambda_{QCD}}$ is a perturbative artifact
First of all, observe that the actual value of $\Lambda_{QCD}$ depends on the order in perturbation theory (and renormalization scheme) in which you are computing the running coupling. By adding higher-order terms to the beta function $\beta(\alpha_{s})$ you change the shape of the running coupling $\alpha_{s}(Q^2)$ which solves the Callan-Symanzik equation for the coupling; therefore the position of the Landau pole is shifted order by order and you have one $\Lambda_{QCD}$ for each non-trivial order in perturbation theory (each of which, except the first and second, depend on the renormalization scheme). You may expect adding orders in perturbation theory to shift the Landau pole to a lower position in momentum space, so as to enable the access to physics at lower and lower momenta. However, this is not what happens: at least to five loops (and for a sufficiently small number of fermions) the higher-order coefficients of the beta function are negative as much as the one-loop order coefficient $\beta_{0}$, so that the derivative of the running coupling is more and more negative and $\alpha_{s}(Q^{2})$ diverges earlier in momentum space ($\Lambda_{QCD}$ is shifted to higher momenta). Of course, there may be an order at which the coefficients change sign and the coupling is allowed to decrease. I must say that currently the evidence is not in favor of this behavior. Setting aside this issue, the point I want to make here is that $\Lambda_{QCD}$ is an intrinsically perturbative scale: it is defined in the context of perturbation theory and has different values at different perturbative orders (and in different renormalization schemes).
The second reason to believe that $\Lambda_{QCD}$ is an artifact of perturbation theory is that QCD is expected to describe the physics of the strong interactions down to zero momentum. Therefore, if the $\alpha_{s}$ that appears in the QCD action is to have any meaning at low momenta, it simply cannot have a Landau pole. Today we know that QCD describes the strong interactions also in the non-perturbative regime thanks to lattice QCD, which was able (for instance) to predict the masses of the meson octet, baryon decuplet and more to an astonishing degree of accuracy. Since lattice QCD exploits an intrinsically non-perturbative approach to computations in QCD, $\Lambda_{QCD}$ is not part of the definition of the theory. Indeed, in lattice QCD a running coupling is not even required to exist. Nonetheless, many definitions can be given of $\alpha_{s}(Q^{2})$ in the lattice framework, all of which must reduce to the standard one in the UV.
$$
$$
An enlightening definition of $\boldsymbol{\alpha_{s}(Q^{2})}$
The definition of $\alpha_{s}(Q^{2})$ which, in my opinion, sheds more light on the interpretation of $\Lambda_{QCD}$ is the one given in the Landau gauge and Taylor scheme, namely
$$
\alpha_{s}(Q^{2})=\alpha_{s}(Q^{2}_{0})\,J(Q^{2};Q^{2}_{0})\,\chi^{2}(Q^{2};Q^{2}_{0}),
$$
where $J(Q^{2};Q^{2}_{0})$ and $\chi(Q^{2};Q^{2}_{0})$ are the gluon and ghost dressing functions renormalized at the scale $Q^{2}_{0}$:
$$
J(Q^{2};Q^{2}_{0})=Q^{2}D(Q^{2};Q^{2}_{0}),\\ \chi(Q^{2};Q^{2}_{0})=Q^{2}G(Q^{2};Q^{2}_{0}),
$$
with $D(Q^{2};Q^{2}_{0})$ and $G(Q^{2};Q^{2}_{0})$ the transverse-gluon and ghost propagators renormalized at the scale $Q^{2}_{0}$ (observe that $J(Q^{2}_{0};Q^{2}_{0})=\chi(Q^{2}_{0};Q^{2}_{0})=1$ by definition). This definition is suitable both in a perturbative and in a non-perturbative setting, since the propagators can be computed in both. In the Landau gauge, it is equivalent to the standard definition of $\alpha_{s}(Q^{2})$ up to two loops. For instance, to one loop and in the Landau gauge, one can compute that
$$
J(Q^{2};Q^{2}_{0})=\left[\frac{\alpha_{s}(Q^{2})}{\alpha_{s}(Q^{2}_{0})}\right]^{13/22}\ ,\qquad \chi(Q^{2};Q^{2}_{0})=\left[\frac{\alpha_{s}(Q^{2})}{\alpha_{s}(Q^{2}_{0})}\right]^{9/44}\ ,
$$
where $\alpha_{s}(Q^{2})$ is the ordinary one-loop running coupling.
On the lattice, one can compute the gluon and ghost propagators and take the product of their dressing functions to obtain a non-perturbative version of $\alpha_{s}(Q^{2})$. The result is contained for instance in Fig. 4 of this article (the computations are made without quarks, but the conclusions are the same). As you can see, on the lattice the Taylor-scheme $\alpha_{s}(Q^{2})$ has no Landau pole: somewhat below 1 GeV (in the absence of quarks), the pole is replaced by a maximum. Moreover, at lower momenta the running coupling decreases until it goes to zero at zero momentum (don't be fooled by this, at zero momentum there can be something else which blows up, giving rise to finite effects!).
This is an example of a running coupling, computed non-perturbatively, which is finite in the IR. What role does $\Lambda_{QCD}$ play in this setting? By itself, none at all. Nonetheless, at high energies the Taylor-scheme coupling computed on the lattice reduces to the standard running coupling. Therefore the high-energy behavior of the Taylor-scheme coupling can indeed be parametrized by the curve (approximating to one loop)
$$
\alpha_{s}(Q^{2})=\frac{4\pi}{\beta_{0}\ln(Q^{2}/\Lambda^{2}_{QCD})}.
$$
Here however $\Lambda_{QCD}$ is a fitting parameter, rather than the position of a pole.
$$
$$
What's going on with the Taylor coupling? Mass effects on the running of the strong coupling
At this point you might be wondering why and how does the Landau pole disappear from the non-perturbative running coupling. In the context of the Taylor scheme and Landau gauge, this question admits a fairly straightforward answer: the finiteness of the coupling can be viewed as being caused by mass effects over its running. In order to illustrate this point I will take QED as an example.
In high-but-not-too-high-energy QED the running coupling can be expressed as
$$
\alpha(Q^{2})=\frac{4\pi}{\beta_{0}\ln(\Lambda^{2}/Q^{2})}\qquad(\beta_{0}>0),
$$
where $\Lambda\sim 10^{286}$ eV can be defined in analogy to $\Lambda_{QCD}$. In the $Q^{2}\to 0$ limit, this expression would imply $\alpha(Q^{2})\to0$, which however is not the correct result. This is because the expression given above does not take into account the mass effects on the running of the coupling due to the electron mass $m_{e}$ being non-zero (recall that most of the elementary derivations of the beta functions go like "Let us suppose that all the masses can be set to zero, then ..."). The correct result instead is
$$
\alpha(Q^{2})=\alpha(Q^{2}_{0})J(Q^{2};Q_{0}^{2})=\frac{\alpha(Q_{0}^{2})}{1-\Pi(Q^{2};Q_{0}^{2})}
$$
where $\Pi(Q^{2};Q_{0}^{2})$ is the photon polarization renormalized at $Q_{0}^{2}$ (notice the similarity of the above equation with the definition of the strong coupling in the Taylor scheme). This expression yields a finite, non-zero result in the limit $\alpha(Q^{2})\to0$, and more generally an IR behavior for $\alpha(Q^{2})$ which is not simply logarithmic.
More generally, at low momenta, one must take into account the mass effects. You might expect that I'm referring to the quarks' masses, as I did above for the electron's mass.  However, I'm not. What I'm talking about is the gluon mass. Indeed, it has now been established that at low energies, due to non-perturbative effects, the gluons acquire a dynamically generated mass. This mass is not expected to explicitly break gauge invariance (although it might be caused by some form of spontaneous symmetry breaking), so it is somewhat a "safe" mass, unlike an explicit mass term in the QCD Lagrangian. At high energies the gluon mass, which is a function of momentum, decreases, until it becomes negligible and the ordinary massless gluons are recovered. The dynamical generation of a mass for the gluons affects the form of the transverse gluon propagator: instead of growing to infinity as $p\to 0$ as would happen for a massless propagator, the gluon propagator saturates to a finite value (see e.g. Fig. 1 in the article I've already cited).
In the context of the Taylor scheme, the existence of a non-perturbative gluon mass scale modifies the form of the beta function with respect to the naive expectations: if there exists an intrinsic mass scale in the theory, then the beta function coefficients are allowed to depend on momentum, rather than being constants. The specific form of these coefficients is framework-dependent, but the general idea is that the gluon mass screens the coupling from becoming infinite by reducing the value of the beta function at small scales: smaller beta implies slower running, hence possibility to avoid the Landau pole.
The results I'm describing cannot be obtained in ordinary perturbation theory: dynamical mass generation for the gluons cannot be described in an ordinary perturbative setting due to perturbative constraints imposed by gauge invariance. Nonetheless, they are currently accepted result which come from lattice studies and other numerical approaches such as those which use the Schwinger-Dyson Equations. Some analytic approaches also managed to obtain similar results.
$$
$$
Conclusions
In non-perturbative formulations of QCD $\Lambda_{QCD}$ does not play a prominent role (if any) in the definition of the running coupling. At best it has the role of a fitting parameter for the high-energy behavior of the coupling. In renormalization schemes such as the Taylor scheme the running coupling can actually be computed in the non-perturbative regime and shown to remain finite. The mass effects caused by the dynamical generation of a mass for the gluons may be responsible for the finiteness of the coupling (this is most certainly true in the Taylor scheme, whereas in other schemes the issue is still open).
$$
$$
Sidenotes


*

*Yes, $T^{\mu}_{\mu}\sim \beta F^{2}$ is valid to all orders.

*Observe that $T^{\mu}_{\mu}$ is RG-invariant and scheme-independent, so $M^{2}$ also is. The product $\beta F^{2}$ is RG-invariant and scheme-independent, unlike the two factors taken separately.

A: I do not pretend that I have the full and right understanding of the problem. The $\Lambda_{QCD}$ is some rough estimate of the energy scale, beyond which physics becomes strongly coupled. Being a perturbative estimate, it can predict the breaking of perturbation series, but solely non-perturbative phenomena, such as instantons, are out of reach. 
As for the trace anomaly, there holds non-renormalization theorem, so this expression is exact, when the exact $\beta$-function is inserted in front of the expression. I've found this reference, maybe you will find it helpful about the problem concerned - https://arxiv.org/abs/1202.1514. It seems, that the results are justified by making an analogy with supersymmetric theories.  
The estimates for $\Lambda_{QCD}$ in literature, therefore, vary significantly : $100-300$ MeV. It is often said, that it is just a some typical scale of order of hadron masses - $\pi, K$ - mesons. 
A: The natural answer to your question is provided by lattice QCD. Like the classical action of QCD (with massless quarks) the lattice (Wilson) action has no dimensionful parameters, only a dimensionless coupling constant $g$. Masses of hadrons are extracted from the exponential decay of correlation functions, so they are expressed in units of the inverse lattice spacing. So a lattice calculation gives the mass of the proton in the form
$$
  M_p={\it const}(g)a^{-1}
$$
which cannot be compared directly to experiment. We compare to experiment by taking the continuum limit $a\to 0$. But $a$ is not a parameter in the dimensionless Wilson action. So what we actually do is to take $g\to 0$, and express $a$ as a function of $g$ using the beta function. This brings $\Lambda_{QCD}$ into the game, because
$$
 a\Lambda_{QCD} = \exp(-1/(2b_0g^2))
$$
(up to higher-loop corrections, which are relevant $a$ is not very small). This fixes the functional form of ${\it const(g)}$, and allows the lattice to determine
$$
 M_p = c_p \Lambda_{QCD}
$$
where $c_p$ is a numerical constant that can be determined by Monte Carlo simulations. This is, of course, what we call dimensional transmutation. We replaced a dimensionless constant, $g$, by a dimensionful one, $\Lambda_{QCD}$.
P.S.: Note that this does not contradict the statement that $\Lambda_{QCD}$ is scheme dependent. I can extract $\Lambda_{QCD}$ from experiment (say the R-ratio) using a particular scheme, for example $\bar{MS}$. I can then do a controlled, perturbative, calculation to relate $\Lambda_{\bar{MS}}$ to $\Lambda_{lat}$, the Lambda parameter for the Wilson lattice regulator. After computing $c_p$ on the lattice I can predict, given the measured R-ratio, the mass of the proton.
P.P.S.: My point in this discussion is to argue that there are no new scales in the infrared, beyond the scale $\Lambda_{QCD}$ that is generated in the UV. In particular, I can take the measured $\alpha_s(M_Z)$ (or the R-ratio, etc) and predict (using the lattice) all non-perturbative quantities, including the string tension, the proton mass, etc.
It is correct that in practice (for computational convenience) most lattice calculations fix the lattice spacing by computing a physical observable, like $f_\pi$ or the string tension. However, they still have to check that the continuum limit is correct. In particular, they need to show that as $a\to 0$ the  coupling goes to zero as dictated by the beta function. This is indeed what Bali and Schilling show. Note that Bali and Schilling also show that I can take the measured $\alpha_s(M_Z)$ and compute the string tension (or the other way around). Finally, Bali and Schilling use the fact that the relation ship between $\Lambda_{\bar{MS}}$ and $\Lambda_{lat}$ is perturbative, and can be determined analytically (as shown by Dashen and Gross).
