Self-energy of a $D$-dimensional statistical mechanical model We consider a $D$ dimensional statistical mechanical model whose partition function is given by 
$$\begin{aligned}
\mathcal{Z} &=\int \mathcal{D} \phi(x) \exp \left(-\mathcal{S}_{\phi}\right) \\
\mathcal{S}_{\phi} &=\int d^{D} x\left\{\frac{1}{2}\left[\left(\nabla_{x} \phi\right)^{2}+r \phi^{2}(x)\right]+\frac{u}{4 !}\left(\phi^{2}(x)\right)^{2}\right\}
\end{aligned}$$
We take a fourier transform of the field and so obtain the action as
$$\begin{aligned}
\mathcal{S}_{\phi}=& \frac{1}{2} \int \frac{d^{D} k}{(2 \pi)^{D}}\left|\phi(k)\right|^{2}\left(k^{2}+r\right) \\
&+\frac{u}{4 !} \int \frac{d^{D} k}{(2 \pi)^{D}} \frac{d^{D} q}{(2 \pi)^{D}} \frac{d^{D} p}{(2 \pi)^{D}} \phi(k) \phi(q) \phi(p) \phi(-k-p-q)
\end{aligned}$$
The susceptibility at $u=0$ is $\chi_0(k)=\frac{1}{k^2+r}$. We perturbatively treat the quartic coupling and so obtain the susceptibility as $$\chi(k)=\frac{1}{1 / \chi_{0}(k)-\Sigma(k)}=\frac{1}{k^{2}+r-\Sigma(k)}$$ where $\Sigma(k)$ is the self energy. Upto a  first order in $u$ we get the self energy as 
$$\Sigma(k)= \frac{-u}{2} \int \frac{d^{D} p}{(2 \pi)^{D}} \frac{1}{p^{2}+r}$$
I don't understand how the self energy and the susceptibility was evaluated. The reference is chapter 3 of Quantum Phase Transitions second edition Subir Sachdev.
 A: Following the comments, it seems that you have understood up to Eq. (3.36) in Sachdev. That is, defining the theory with the partition function
$$
\mathcal{Z} = \int \left( \prod_{i}  \frac{d y_i}{\sqrt{2 \pi}} \right) \exp \left( - \frac{1}{2} \sum_{ij} y_i A_{ij} y_j - \frac{u}{4!} \sum_i y_i^4 \right),
$$
the two-point function $C_{ij} = \langle y_i y_j \rangle$ is equal to
$$
C_{ij} = A^{-1}_{ij} - \frac{u}{2} \sum_{k} A_{ik}^{-1} A^{-1}_{kk} A^{-1}_{kj} + \cdots
$$
to leading order in $u$. (Sachdev gives expressions to order $u^2$, but we won't need this to derive the term you have asked about.)
Now we want to calculate $C^{-1}_{ij}$ perturbatively in $u$. First of all, since $C_{ij} = A^{-1}_{ij}$ when $u=0$, clearly we can write
$$
C^{-1}_{ij} = A_{ij} + u \, \kappa_{ij} + \cdots
$$ 
where we are dropping terms of order $u^2$ or higher, and $\kappa_{ij}$ is a $u$-independent matrix to be determined. Now by definition of the inverse,
$$
C_{ij} C^{-1}_{jm} = \delta_{im} .
$$
We now just plug in the above expressions for $C_{ij}$ and $C^{-1}_{ij}$, finding
$$
C_{ij} C^{-1}_{jm} = \left( A^{-1}_{ij} - \frac{u}{2} \sum_{k} A_{ik}^{-1} A^{-1}_{kk} A^{-1}_{kj} \right) \left( A_{jm} + u \, \kappa_{jm} \right) + \cdots \\ = \delta_{im} + u A^{-1}_{ij} \kappa_{jm} - \frac{u}{2} \sum_k A^{-1}_{ik} A^{-1}_{k k} A^{-1}_{kj} A_{jm} + \cdots
$$
Once again, I'm always ignoring terms higher-order in $u$. We set this equal to $\delta_{i m}$, and use the fact that $A^{-1}_{kj} A_{jm} = \delta_{km}$ by the definition of inverses. At this point the algebra is especially simple, and you should find
$$
\kappa_{ij} = \frac{1}{2} \delta_{ij} A^{-1}_{kk}.
$$
Alternatively, using the notation of (3.37),
$$
C^{-1}_{ij} = A_{ij} - \Sigma_{ij}
$$
and comparing with our previous expressions,
$$
\Sigma_{ij} = - \frac{u}{2} \delta_{ij} A^{-1}_{kk}.
$$
At this point, we need to determine what $A^{-1}_{ij}$ is for the continuum model you have written down. There are a few approaches here. One is to discretize the model, in which case everything looks just as it does above, and one gets your final expression after taking the continuum limit of the final discrete expression. Alternatively, there are effectively "continuum" expressions for the above equations, where you replace
$$
A_{ij} A^{-1}_{jm} = \delta_{im}
$$
with
$$
\int d^D y \, K(x,y) G(y,z) = \delta^D(x - z),
$$
where we have the Dirac delta on the right-hand side, $G(x,y)$ is called the Green's function, and and $K(x,y)$ is the function appearing in the continuum action,
$$
\mathcal{S}_0 = \frac{1}{2} \int d^D x \, d^D y \, \phi(x) K(x,y) \phi(y).
$$
For your initial action, we have $K(x,y) = \delta^{D}(x - y) \left( - \nabla^2_x + r \right)$, and the Green's function satisfies
$$
\left( - \nabla^2_x + r \right) G(x,y) = \delta^D(x - y).
$$
This is solved using a Fourier transformation, and we get
$$
G(x,y) = \int \frac{d^D p}{(2 \pi)^D} \frac{e^{i p \cdot (x - y)}}{p^2 + r}
$$
Now, using this expression in the continuum limit of the above equation for $\Sigma$, we find
$$
\Sigma(x,y) = - \frac{u}{2} \delta^D(x - y) \int \frac{d^D p}{(2 \pi)^D} \frac{1}{p^2 + r}.
$$
Or, after a Fourier transformation,
$$
\Sigma(k) = - \frac{u}{2}  \int \frac{d^D p}{(2 \pi)^D} \frac{1}{p^2 + r}.
$$
