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.
