Let us start with time-independent Brillouin-Wigner (BW) perturbation theory for a non-degenerate isolated eigenstate, i.e. it is assumed to not be part of a continuous spectrum. Here we will consider a (not necessarily convergent) formal power series in a coupling constant $\lambda$.
Using standard notation we have
$$\begin{array}{lrcll}\text{Hamiltonian}& \hat{H}&=&\hat{H}_0+\lambda\hat{V}&\cr
\text{Unperturbed eigenstate}& \hat{H}_0|\psi_0\rangle&=& E_0|\psi_0\rangle,& \langle\psi_0|\psi_0\rangle~=~1,\cr
\text{Perturbed eigenstate}& \hat{H}|\psi\rangle&=& E|\psi\rangle,&\cr
\text{Energy shift}&\Delta&:=&E-E_0&\cr
& &=&\sum_{n\in\mathbb{N}} E_n\lambda^n&~=~{\cal O}(\lambda),\cr
&E&=&\sum_{n\in\mathbb{N}_0} E_n\lambda^n,& \cr
\text{Projection operator}& \hat{P}&:=&|\psi_0\rangle\langle\psi_0|,& [\hat{P},\hat{H}_0]~=~0,\cr
\text{Projection operator}& \hat{Q}&:=&\hat{\bf 1}-\hat{P}~=~&(E-\hat{H}_0)\hat{R},\cr
\text{Resolvent}& \hat{R}&:=&\frac{1}{E-\hat{H}_0}\hat{Q}~=~&\frac{1}{E_0-\hat{H}_0+\Delta}\hat{Q}\cr
& &=&\frac{E_0-\hat{H}_0}{E_0-\hat{H}_0+\Delta}\hat{R}_0
~=~&\frac{1}{1+\hat{R}_0\Delta}\hat{R}_0\cr
& &=&\sum_{m\in\mathbb{N}_0}(-\Delta)^m\hat{R}_0^{m+1},\cr
\text{Resolvent}& \hat{R}_0&:=&\frac{1}{E_0-\hat{H}_0}\hat{Q}.
\end{array}\tag{1}$$
We now assume that the (not necessarily normalized) perturbed state $|\psi\rangle$ has a non-zero overlap
$$ \langle\psi_0|\psi\rangle~\neq~0 \tag{2}$$
with the unperturbed state $|\psi_0\rangle$. We may then rescale $|\psi\rangle$ with a non-zero (possibly $\lambda$-dependent) constant:
$$ \hat{P}|\psi\rangle~=~|\psi_0\rangle
\quad\Rightarrow\quad
\langle\psi_0|\psi\rangle~=~1
\quad\wedge\quad
\langle\psi|\psi\rangle~>~1. \tag{3}$$
In this way we obtain a so-called intermediate normalization of the overlap (3).
Then
$$|\psi\rangle-|\psi_0\rangle~\stackrel{(3)}{=}~\hat{Q}|\psi\rangle
~\stackrel{(1)}{=}~\hat{R}(E-\hat{H}_0)|\psi\rangle
~\stackrel{(1)}{=}~\hat{R}\lambda\hat{V}|\psi\rangle,\tag{4}$$
which leads to a fixed-point equation for the perturbed state
$$ |\psi\rangle~\stackrel{(4)}{=}~|\psi_0\rangle+\hat{R}\lambda\hat{V}|\psi\rangle~\stackrel{(4)}{=}~\sum_{n\in\mathbb{N}_0}(\hat{R}\lambda\hat{V})^n|\psi_0\rangle~=~\frac{1}{1-\hat{R}\lambda\hat{V}}|\psi_0\rangle.\tag{5}$$
This leads to the BW fixed-point equation for the shifted energy
$$\begin{align}\Delta~\stackrel{(1)}{:=}~&E-E_0\cr
~\stackrel{(3)}{=}~&\langle\psi_0|(E-E_0)|\psi\rangle\cr
~\stackrel{(1)}{=}~&\langle\psi_0|(E-\hat{H}_0)|\psi\rangle\cr
~\stackrel{(1)}{=}~&\langle\psi_0|\lambda\hat{V}|\psi\rangle\cr
~\stackrel{(5)}{=}~&\langle\psi_0|\lambda\hat{V}\frac{1}{1-\hat{R}\lambda\hat{V}}|\psi_0\rangle\cr
~\stackrel{(5)}{=}~&\sum_{n\in\mathbb{N}_0}\langle\psi_0|\lambda\hat{V}(\hat{R}\lambda\hat{V})^n|\psi_0\rangle\cr
~\stackrel{(1)}{=}~&\langle\psi_0|\lambda\hat{V}\frac{1}{1-\frac{1}{1+\hat{R}_0\Delta}\hat{R}_0\lambda\hat{V}}|\psi_0\rangle\cr
~\stackrel{(1)}{=}~&\langle\psi_0|\lambda\hat{V}\frac{1}{1-\sum_{m\in\mathbb{N}_0}(-\Delta)^m\hat{R}_0^{m+1}\lambda\hat{V}}|\psi_0\rangle\cr
~=~&\sum_{n\in\mathbb{N}_0}\lambda^{n+1}\sum_{m_1,\ldots,m_n\in\mathbb{N}_0}(-\Delta)^{\sum_{i=1}^nm_i}\cr
&\underbrace{\langle\psi_0|\hat{V}\hat{R}_0^{m_1+1}\hat{V}\hat{R}_0^{m_2+1}\hat{V}\ldots\hat{V}\hat{R}_0^{m_n+1}\hat{V}|\psi_0\rangle}_{~=:~ \langle m_1,m_2,\ldots,m_n\rangle}\cr
~\stackrel{(1)}{=}~&\sum_{n\in\mathbb{N}_0}\lambda^{n+1}
\sum_{M\in\mathbb{N}_0}(-1)^M
\sum_{m_1,\ldots,m_n\in\mathbb{N}_0}^{\sum_{i=1}^nm_i=M}\left(\sum_{n_0\in\mathbb{N}} E_{n_0}\lambda^{n_0}\right)^M \langle m_1,m_2,\ldots,m_n\rangle.
\end{align}\tag{6}$$
We may then systematically obtain a formula for the $N$th-order energy correction
$$\begin{align}E_N~\stackrel{(6)+(8)}{=}~&
\sum_{n=0}^{N-1}\sum_{M=0}^{N-(n+1)}(-1)^MP_{N-(n+1),M}\cr
&\sum_{m_1,\ldots,m_n\in\mathbb{N}_0}^{\sum_{i=1}^nm_i=M}
\underbrace{\langle\psi_0|\hat{V}\hat{R}_0^{m_1+1}\hat{V}\hat{R}_0^{m_2+1}\hat{V}\ldots\hat{V}\hat{R}_0^{m_n+1}\hat{V}|\psi_0\rangle}_{~=:~ \langle m_1,m_2,\ldots,m_n\rangle}.
\end{align}\tag{7}$$
in Rayleigh-Schrödinger (RS) perturbation theory.
Here the polynomials
$$P_{N,M}(E_1,\ldots,E_{N-M+1})~:=~\sum_{n_1,\ldots,n_M\in\mathbb{N}}^{\sum_{i=1}^M n_i=N}\prod_{j=1}^ME_{n_j}\tag{8}$$
have a generating function
$$ P_{N,M}\lambda^Nq^M~=~\sum_{M\in\mathbb{N}_0}(q\Delta)^M~=~\frac{1}{1-q\Delta}.\tag{9}$$
The first few polynomials read
$$\begin{align}
P_{N,0}~=~&\delta_N^0,\cr
P_{N,1}~=~&E_N,\cr
P_{N,2}~=~&\sum_{n-1}^{N-1}E_nE_{N-n},\cr
~\vdots&\cr
P_{N,N}~=~&E_1^N.
\end{align}\tag{10}$$
One can check that the polynomial in eq. (7) only depends on the previous energy corrections $E_1,\ldots, E_{N-2}$, so that eq. (7) constitutes a recurrence relation.
The first few orders of energy corrections read
$$\begin{align}
E_1~\stackrel{(7)}{=}~&P_{0,0}\langle\rangle
~=~ \langle\psi_0|\hat{V}|\psi_0\rangle, \cr
E_2~\stackrel{(7)}{=}~&P_{0,0}\langle 0\rangle
~=~\langle\psi_0|\hat{V}\hat{R}_0\hat{V}|\psi_0\rangle,\cr
E_3~\stackrel{(7)}{=}~&P_{0,0}\langle 00\rangle-P_{1,1}\langle 1\rangle\cr
~=~&\langle\psi_0|\hat{V}\hat{R}_0\hat{V}\hat{R}_0\hat{V}|\psi_0\rangle-E_1\langle\psi_0|\hat{V}\hat{R}_0^2\hat{V}|\psi_0\rangle,\cr
E_4~\stackrel{(7)}{=}~&P_{0,0}\langle 000\rangle-2P_{1,1}\langle 10\rangle -P_{2,1}\langle 1\rangle +P_{2,2}\langle 2\rangle\cr
~=~&\langle\psi_0|\hat{V}\hat{R}_0\hat{V}\hat{R}_0\hat{V}\hat{R}_0\hat{V}|\psi_0\rangle
-2E_1\langle\psi_0|\hat{V}\hat{R}_0^2\hat{V}\hat{R}_0\hat{V}|\psi_0\rangle\cr
&-E_2\langle\psi_0|\hat{V}\hat{R}_0^2\hat{V}|\psi_0\rangle
+E_1^2\langle\psi_0|\hat{V}\hat{R}_0^3\hat{V}|\psi_0\rangle,\cr
E_5~\stackrel{(7)}{=}~&P_{0,0}\langle 0000\rangle -P_{1,1}(2\langle 100\rangle+\langle 010\rangle) -2P_{2,1}\langle 10\rangle\cr
&+P_{2,2}(2\langle 20\rangle+\langle 11\rangle) -P_{3,1}\langle 1\rangle +P_{3,2}\langle 2\rangle-P_{3,3}\langle 3\rangle,\cr
E_6~\stackrel{(7)}{=}~&P_{0,0}\langle 00000\rangle -P_{1,1}(2\langle 1000\rangle+\langle 0100\rangle)\cr
&-P_{2,1}(2\langle 100\rangle+\langle 010\rangle)\cr
&+P_{2,2}(2\langle 200\rangle+\langle 020\rangle+2\langle 110\rangle+\langle 101\rangle)\cr
&-2P_{3,1}\langle 10\rangle +P_{3,2}(2\langle 20\rangle+\langle 11\rangle) -2P_{3,3}(\langle 30\rangle+\langle 21\rangle)\cr
&-P_{4,1}\langle 1\rangle +\underbrace{P_{4,2}}_{=2E_3E_1+E_2^2}\langle 2\rangle
-\underbrace{P_{4,3}}_{=3E_2E_1^2}\langle 3\rangle +P_{4,4}\langle 4\rangle,
\end{align}\tag{11}$$
and so forth.
