To the best of my knowledge of the literature on this topic, the answer is: not really (a few exceptions appear below). Let me first give a rigorous statement of the problem: the physical equation $\beta = 0$ can be expressed as
$\sum_{k=0}^{\infty} \epsilon^k \beta_k = 0$,
where $\epsilon$ is a physical parameter one imagines is "small," and where $\beta_k$
is a symmetric two tensor depending on a Riemannian metric which arises as a certain universal expression in the curvature tensor and its derivatives, and which moreover obeys a natural scaling law and a natural "homogeneity" in terms of the number of derivatives of $g$ which are required to express it. These terms $\beta_k$ are computable in theory, although my understanding is that this has only been done on a piecemeal basis for the first few terms. There are many references one can easily find which do this for this and other sigma models. The first two terms are $\beta_0 = \mbox{Vol}(M) g$, $\beta_1 = \mbox{Rc}$, the term $\beta_2$ is a quadratic expression in the curvature tensor. Considering the equation up to order $k$ in the $\epsilon$ power expansion is referred to in physics literature as "up to $k$-loops."
So, with this background, there are at least two interesting questions one can ask:
1) Given $N \in \mathbb N \cup \{\infty\}$, can one construct on a given manifold a solution to $\sum_{k=0}^N \epsilon^k \beta_k = 0$?
2) Can one construct on a given manifold a one-parameter family of Riemannian metrics satisfying $\frac{\partial g}{\partial t} + \sum_{k=0}^N \epsilon^k \beta_k = 0$?
The original question was related to 1) above, but 2) is as important in the physics literature (AFAIK). Of course, the case $N=1$ of question 1) corresponds to solving the usual Riemannian Einstein equation (although one may need to justify ignoring the term $k=0$ in some cases, which physicists have ways of doing). In the Kähler setting the metric can come from the Calabi-Yau theorem as the poster stated. For question 2), the case $N=1$ corresponds to the (properly normalized) Ricci flow equation.
Now that I have properly stated the question, let me just say that I focus on studying elliptic/parabolic equations on manifolds, and for a while got interested in exactly this question. After much literature digging, I found only a few rigorous mathematical results:
http://arxiv.org/abs/0904.1241
http://arxiv.org/abs/1108.0526
http://arxiv.org/abs/1205.6507
The first paper considers the case $N=2$ of question 2). In this case the operator is still second-order, and so in special settings the equation can be rendered parabolic, and then known techniques can be applied. This paper is certainly interesting, although the techniques cannot really be extended to $N > 2$. The last two papers address the homogeneous setting.
Moving beyond $N = 2$ is, IMHO, extremely difficult from a PDE point of view, since one has very little information on the form of the terms $\beta_k$ beyond the general qualitative statements I made above. For instance, an interesting side question is: do there exist arbitrarily large values of $k$ for which $\beta_k$ corresponds to an elliptic operator of a Riemannian metric? This is true for $k=1$, but I believe it is very unlikely to be true for higher $k$ since I believe it is "known" in the physics sense that all of the terms arising in $\beta_k$ for $k > 1$ are at least quadratic combinations of curvature and derivatives. I.e., one cannot expect a term like $\Delta \mbox{Rc}$ to show up, which would be elliptic.
One can imagine considering the case $N > 2$ on homogeneous spaces, where in principle the equations reduce to a system of ODE, but again one does not have much in the way of understanding the form of the general term $\beta_k$.
It is also not inconceivable that if one starts with say a Ricci-flat metric, and is allowed to choose $\epsilon$ very small with respect to this metric, that one can perturb it to a solution to the higher order equations, but again without at least a little more hard data on the form of the $k$-th term this seems tricky.
EDIT: In response to Robert's comments, yes, I meant $\mbox{Vol}(M) g$ (fixed above). Also, let me give a reference:
http://www.sciencedirect.com/science/article/pii/0550321389904227
This paper computes the terms $\beta_3$ and $\beta_4$ (The formula for $\beta_4$ covers half a page of terms, and each of these is shorthand for a complicated curvature expression already. Written out fully would take 2 pages at least!). As far as whether one can expect the $\beta$ functions to be simpler on a Kähler manifold, I can't speak to the physics well enough to really say. From a mathematical point of view, I don't really see any advantage beyond the fact that of course each term could be expressed in terms of a Kähler potential.
The case of symmetric spaces is probably more tractable. Looking for solutions within the class $\nabla \mbox{Rm} \equiv 0$ could certainly simplify the terms $\beta_k$ greatly. My understanding of how these terms are derived is extremely fuzzy, but let me speak anyways: I believe the different summands in $\beta_k$ correspond in some sense to different sized "loop diagrams" (see the above reference for some examples), and moreover the qualitative behavior of such terms (i.e. number of derivatives of $\mbox{Rm}$ which appear is (somehow) related to the topology of the diagram. One main difficulty in performing these calculations to derive the terms $\beta_k$ is the sheer number of diagrams which must be considered, which grows wildly with $k$. If, on the other hand, one could a priori rule out a large number of such diagrams (by assuming $\nabla \mbox{Rm} \equiv 0$), perhaps these calculations become more tractable, or at least some stronger qualitative statements could be made.
Lastly, the quantity $\epsilon$ (which, also note is usually $\alpha'$ in physics literature) is, I believe, meant to be small but fixed. Physically it apparently is meant to represent the "string tension."
Best Answer
The reason for this is Allard's regularity theorem.
Roughly speaking Allard's theorem says that if near a point of the support of a stationary varifold the varifold has unit density and area close to that of the ball of the appropriate dimension (for a 2-varifold it would be area of a disk) then the support of the varifold is smooth at that point.
More precisely, there is an $\epsilon>0$ and $r_0>0$ (depending on the ambient geometry and dimension of the varifold). So that if $\Sigma$ is an stationary $m$-varifold in $M$ , a Riemannian manifold, and a point $p\in spt \Sigma$ satisfies $$\mathcal{H}^m(B_r(p)\cap \Sigma) \leq \omega_m r^m(1+\epsilon)$$ for $r\leq r_0 $ then $spt \Sigma$ is smooth near $p$. Here $B_r(p)$ is the $r$-ball in $M$ and $\omega_m$ is the volume of the unit ball in $\mathbb{R}^m$.
Allard's proof is unfortunately quite technical -- a good reference is Leon Simon's (sadly) hard to find book "Lectures on Geometric Measure Theory". Luckily, for your purpose there is a simpler version with a very easy proof due to Brian White (see here for the paper).
Specifically, suppose you have instead of being a stationary varifold you know that $\Sigma$ is a smooth minimal surface. If $p$ is a point of $\Sigma$ so that $$\mathcal{H}^m(B_r(p)\cap \Sigma) \leq \omega_m r^m(1+\epsilon)$$ then one has $$|A|(p)\leq r^{-1}$$ here $|A|$ is the norm of the second fundamental form. In other words you obtain a quantitive bound on curvature.
How does this relate to your question? Well you have that $\Sigma_k$ converge to $\Sigma$ as Radon measures. Let $p\in \Sigma$ this convergence implies that $$\mathcal{H}^m(\Sigma_k\cap B_r(p))\to \mathcal{H}^m(\Sigma\cap B_r(p))$$ (it is worth noting that we are using that the surfaces are area minimizing to ensure the convergence is with multiplicity one). Now since $\Sigma$ is smooth near $p$ it is locally modelled on a flat plane. In other words, there is a scale $r$ so that $$\mathcal{H}^m(B_r(p)\cap \Sigma) \leq \omega_m r^m(1+\frac{1}{2}\epsilon)$$ hence by the convergence $$\mathcal{H}^m(B_r(p)\cap \Sigma_k) \leq \omega_m r^m(1+\epsilon)$$ and so since the $\Sigma_k$ are smooth (either a priori or by Allard's full theorem) and the point $p$ was not important we deduce that for some $\delta>0$ we have $$\sup_{B_\delta(p)\cap\Sigma_k} |A|\leq r^{-1}.$$ In other words there is a uniform curvature bound on the $\Sigma_k$. The result then follows from "Standard elliptic PDE" and the Arzela-Ascoli theorem.