Calculating exactly the divergent part of amplitudes at all loop order

Question

Suppose we have the $L$-loop amplitude of the form 

$$\mathcal{I}_L=\int \prod_{i=1}^L \frac{d^D q_i}{(2 \pi)^D} \frac{1}{q_i^2} \frac{1}{(p-\sum_{i=1}^L q_i)^2}.$$

Introducing Feynman parameters to merge the denominators as usual, we may write the amplitude in the form 

$$\mathcal{I}_L = (N-1)! \int_0^1 \prod_{j=1}^{L+1} dx_j \delta \left( \sum_{i=1}^{L+1} x_i-1\right) \int \prod_{i=1}^L \frac{d^D q_i}{(2 \pi)^D} \left[ q_i q_j M_{ij}-2 q_j K_j+J \right]^{-(L+1)}$$

where $M$ is a symmetric matrix. We can evaluate the $q$ integrals using the formula 

$$\int \frac{d^D l}{(2 \pi)^D}\frac{1}{(l^2-\Delta)^N} = \frac{(-1)^Ni}{(4 \pi)^{D/2}}\frac{\Gamma(N-D/2)}{\Gamma(N)}\left( \frac{1}{ \Delta}\right)^{N-D/2},$$

It can be shown that we can write $\mathcal{I}_L$ as 

$$\mathcal{I}_L=\frac{(-1)^L \Gamma(L+1-D)}{(4 \pi)^D} \int_0^1 \prod_{i=1}^{L+1} dx_i \delta \left( \sum_{i=1}^{L+1} x_i-1\right) \frac{\mathcal{U}^{L+1-3D/2}}{\mathcal{F}^{L+1-D}},$$

where we have defined  $\mathcal{U}=\det M$ and $\mathcal{F}=\det M \left(  K_i M^{-1}_{ij} K_j -J \right).$

By setting $D=4-2 \epsilon$, my goal is to calculate the divergent part 

$$\mathcal{I}_L = \frac{c_L}{\epsilon} + \mathcal{O}(\epsilon^0)$$

for any $L$. With that aim in mind, let's first look at $L=2$. Doing the above computations and calculating $\mathcal{U}$ and $\mathcal{F}$ explicitly, amounts to 

$$c_L \propto \int dx_1 dx_2 dx_3 \delta \left( x_1+x_2+x_3-1\right) \frac{x_1 x_2 x_3}{(x_1 x_2 +x_2 x_3 + x_1 x_3)^3}.$$

The form for larger $L$ follows a similar pattern (product of $x_i$ in the numerator to some power, and sum of products of $x_i$ with one missing in each term, just like above), so evaluating the integral for $L=2$ would probably lead the way for a more general result. However, how does one treat integrals such as this? Any ideas about how one might try to evaluate it completely?

Comments

The sweeping under the carpet of divergencies is one of the reasons I never delved deeply into QFT, and therefore I restrict myself to a comment instead of an answer.

The $\delta$-function as such is pernicious evil, obscuring often what is meant. If the above is indeed a volume integral, then the result is 0 (taking the $\delta$-function as what it is usually described to be, infinity at one point and zero everywhere else, and using Lebesgue-integrals). If you take it as a distribution then things should be expressed differently, or one could employ a limiting procedure (which would point a way to a solution, as it makes clearer what is meant).

What probably is meant, however, is that there is a contribution to the integral only if $x_1+x_2+x_3-1=0$, which then suggests that the integral is not intended as a volume integral but as a surface integral. Why not substitute $x_3=1-x_1-x_2$ and integrate over $x_1$and $x_2$ only?
 

---

Thanks for the comment! Indeed, I mean the latter. Of course, substituting for $x_3$ is the first thing I tried, but it is still not very tractable. I left it as is because it has a symmetric form that is similar to a more general expression for all $L$. But yes, it is essentially a surface integral.