Higher orders in perturbation theory I would like to compute an energy level up to many orders in perturbation theory. My difficulty right now is not in the calculation itself but in understanding the algebraic structure of the higher orders (above the second order).
The higher orders contain many terms, and I can't see much rhyme or reason to how they work. Say, according to Wikipedia the fourth order correction is given by 
$$E_n^{(4)} = \frac{V_{nk_4}V_{k_4k_3}V_{k_3k_2}V_{k_2n}}{E_{nk_2}E_{nk_3}E_{nk_4}}-\frac{|V_{nk_4}|^2}{E_{nk_4}^2}\frac{|V_{nk_2}|^2}{E_{nk_2}}-V_{nn}\frac{V_{nk_4}V_{k_4k_3}V_{k_3n}}{E_{nk_3}^2E_{nk_4}}-V_{nn}\frac{V_{nk_4}V_{k_4k_2}V_{k_2n}}{E_{nk_2}E_{nk_4}^2}+V_{nn}^2\frac{|V_{nk_4}|^2}{E_{nk_4}^3}$$ 
with all $k$ indices being summed over so that denominators do not vanish. Fifth order is a nightmare. 
Can these higher order expressions be written in a more systematic way? As sums over integer partitions or set partitions or some other combinatorial structure?
 A: OP asks about the algebraic structure (rather than the actual value) of the $m$'th order $E^{(m)}_n$ of the $n$'th energy level in non-degenerate perturbation theory, see Wikipedia to fix notation.
It is natural to introduce a type of "Feynman diagrams" to indicate the algebraic structure. The actual value is encoded in an integer coefficient/weight in front of the Feynman diagram, which we do not discuss here.$^a$
Feynman rules:


*

*External sources $x$ with same label $n$.

*Vertex $o$ with summation label $k\neq n$, which we later should sum over.

*Horizontal oriented line corresponding to matrix element $V_{k_2k_1}$.

*Non-horizontal oriented line corresponding to $1/E_{nk_1}$ from a vertex $o$ to an external source $x$.

*A vertex $o$ has one incoming and one outgoing horizontal leg and at least one outgoing non-horizontal leg.

*At $m$'th order there is $m$ horizontal lines; $m\!-\!1$ non-horizontal lines, and at most $m\!-\!1$ vertices $o$. 
Example:
  x---<---o---<---o---<---x    x---<---x    x---<---x
         /|\      |
        / | \     |
       x  x  x    x

$\uparrow$ Fig. 1. A typical non-connected Feynman diagram with $m=5$ corresponding to the term
$$V_{nk_2}\frac{1}{E_{nk_2}} \frac{1}{E_{nk_2}}\frac{1}{E_{nk_2}}V_{k_2k_1}\frac{1}{E_{nk_1}}V_{k_1n}~ V_{nn}~V_{nn}.$$
--
$^a$ The sign of the integer coefficient/weight is equal to $-(-1)^{\#(\text{conn. comp.})}$.
A: In case the potential is a polynomial, the perturbative expansion for the energy eigenvalues can be generated using very simple recursion relations as shown first by Bender and Wu in Phys. Rev. 184, 1231 (1969) and Phys. Rev. D7, 1620 (1973), see page 29 and further here for the details. This has allowed people to obtain thousands of terms to test resummation methods for divergent perturbation series. 
