Preliminary
First of all, we work with the mostly-plus metric convention $g^{\mu\nu} = \mathrm{diag}(-1,1,1,1)$. Also, we adopt the dimensional regularization of divergence with taking the spacetime dimension $d=4-2\epsilon$.Dirac Algebra
We first summarize some Dirac algebra in $4$ and $4-2\epsilon$ dimesnion. Followings are contraction formulas in $4$ [$4-2\epsilon$] dimension: \begin{align} \gamma^\mu \gamma_\mu &= 4\ [d], \\ \gamma^\mu \gamma^\nu \gamma_\mu &= -2\gamma^\nu\ [(2-d)\gamma^\nu], \\ \gamma^\mu \gamma^\nu \gamma^\rho \gamma_\mu &= 4g^{\nu\rho}\ [4g^{\nu\rho} - (4-d) \gamma^\nu \gamma^\rho],\\ \gamma^\mu \gamma^\nu \gamma^\rho \gamma^\sigma \gamma_\mu &= -2\gamma^\sigma \gamma^\rho \gamma^\nu\ [-2\gamma^\sigma \gamma^\rho \gamma^\nu + (4-d) \gamma^\nu \gamma^\rho \gamma^\sigma]. \end{align} On the other hand, the trace formulas are independent of the choice of the spacetime dimension and \begin{align} {\rm Tr}[ \gamma^\mu \gamma^\nu ] &= 4 g^{\mu\nu},\\ {\rm Tr}[ \gamma^\mu \gamma^\nu \gamma^\rho \gamma^\sigma ] &= 4(g^{\mu\nu} g^{\rho\sigma} - g^{\mu\rho} g^{\nu\sigma} + g^{\mu\sigma} g^{\nu\rho}),\\ {\rm Tr}[ \gamma^\mu \gamma^\nu \gamma^\rho \gamma^\sigma \gamma_5 ] &= -4i\epsilon^{\mu\nu\rho\sigma}. \end{align}Feynman Parameters
Following equalities using Feynman parameters are useful to perform the $d$-dimensional integral \begin{align} \frac{1}{AB} &= \int_0^1 dxdy\ \delta(1-x-y) \frac{1}{[xA + yB]^2},\\ \frac{1}{AB^n} &= \int_0^1 dxdy\ \delta(1-x-y) \frac{ny^{n-1}}{[xA+yB]^{n-1}}, \\ \frac{1}{A_1\cdots A_n} &= \int_0^1 dx_1 \cdots dx_n\ \delta\left(1-\sum_{i=1}^{n} dx_i\right) \frac{(n-1)!}{\left[\sum_{i=1}^{n} x_i A_i\right]^n}, \end{align}Loop Integral Formulas
Considering the $d$-dimensional integral, we obtain the general formula \begin{align} \mu^{4-d} \int \frac{d^d k}{(2\pi)^d} \frac{k^{2a}}{(k^2 - \Delta)^b} = i(-1)^{a-b} \frac{\mu^{4-d}}{(4\pi)^{d/2}} \frac{1}{\Delta^{b-a-\frac{d}{2}}} \frac{\Gamma (a+\frac{d}{2}) \Gamma (b-a-\frac{d}{2})}{\Gamma (b) \Gamma (\frac{d}{2})}. \end{align} In particular, following examples are convergent when we take $d=4$, \begin{align} \int \frac{d^4 k}{(2\pi)^4} \frac{1}{(k^2 - \Delta)^3} &= \frac{-i}{32\pi^2 \Delta},\\ \int \frac{d^4 k}{(2\pi)^4} \frac{k^2}{(k^2 - \Delta)^4} &= \frac{-i}{48\pi^2} \frac{1}{\Delta}, \\ \int \frac{d^4 k}{(2\pi)^4} \frac{1}{(k^2 - \Delta)^r} &= i\frac{(-1)^r}{(4\pi)^2} \frac{1}{(r-1)(r-2)} \frac{1}{\Delta^{r-2}},\\ \int \frac{d^4 k}{(2\pi)^4} \frac{k^2}{(k^2 - \Delta)^{r+1}} &= i\frac{(-1)^r}{(4\pi)^2} \frac{2}{r(r-1)(r-2)} \frac{1}{\Delta^{r-2}}, \end{align} with $r \geq 3$, while the divergent examples are \begin{align} \mu^{4-d} \int \frac{d^d k}{(2\pi)^d} \frac{1}{k^2 - \Delta} &= \frac{-i \mu^{4-d}}{(4\pi)^{d/2}} \frac{1}{\Delta^{1-d/2}} \Gamma\left( \frac{2-d}{2} \right) \simeq \frac{i\Delta}{16\pi^2} \left[ \frac{1}{\epsilon} + 1 + \ln \frac{\tilde{\mu}^2}{\Delta} \right], \\ \mu^{4-d} \int \frac{d^d k}{(2\pi)^d} \frac{1}{[k^2 - \Delta]^2} &= \frac{i \mu^{4-d}}{(4\pi)^{d/2}} \frac{1}{\Delta^{2-d/2}} \Gamma\left( \frac{4-d}{2} \right) \simeq \frac{i}{16\pi^2} \left[ \frac{1}{\epsilon} + \ln \frac{\tilde{\mu}^2}{\Delta} \right], \\ \mu^{4-d} \int \frac{d^d k}{(2\pi)^d} \frac{k^2}{[k^2 - \Delta]^2} &= -\frac{d}{2} \frac{i \mu^{4-d}}{(4\pi)^{d/2}} \frac{1}{\Delta^{1-d/2}} \Gamma\left( \frac{2-d}{2} \right) \simeq \frac{i\Delta}{16\pi^2} \left[ \frac{2}{\epsilon} + 1 + 2\ln \frac{\tilde{\mu}^2}{\Delta} \right],\\ \mu^{4-d} \int \frac{d^d k}{(2\pi)^d} \frac{k^4}{[k^2 - \Delta]^2} &= \frac{d}{2} \left( 1 + \frac{d}{2} \right) \frac{i\mu^{4-d}}{(4\pi)^{d/2}} \frac{1}{\Delta^{-d/2}} \Gamma \left( -\frac{d}{2} \right) \simeq \frac{i\Delta^2}{16\pi^2} \left[ \frac{3}{\epsilon} + 2 + 3\ln \frac{\tilde{\mu}^2}{\Delta} \right],\\ \mu^{4-d} \int \frac{d^d k}{(2\pi)^d} \frac{k^2}{[k^2 - \Delta]^3} &= \frac{d}{4} \frac{i \mu^{4-d}}{(4\pi)^{d/2}} \frac{1}{\Delta^{2-d/2}} \Gamma\left( \frac{4-d}{2} \right) \simeq \frac{i}{16\pi^2} \left[ \frac{1}{\epsilon} - \frac{1}{2} + \ln \frac{\tilde{\mu}^2}{\Delta} \right], \end{align} with $\tilde{\mu}^2 \equiv 4\pi e^{-\gamma_E} \mu^2$.PaVe Functions
It is often useful to consider the Passarino-Veltoman basis of the one-loop integrals.1 In our case, there are only $4$ linearly independent basis of integrals, i.e., $1$-, $2$-, $3$-, and $4$-point scalar integrals. In the following, we will see them one by one and also see relevant reduction formulas of vector and tensor integrals.1-point function
The only important integral is \begin{align} A(m) \equiv \frac{\mu^{4-d}}{i\pi^2} \int d^dk \frac{1}{k^2-m^2} = m^2 \left[ \frac{1}{\epsilon} + 1 + \ln \frac{\tilde{\mu}^2}{m^2} \right]. \end{align}2-point functions
We define \begin{align} B_0; B_\mu; B_{\mu\nu}(p, m_1^2, m_2^2) \equiv \frac{\mu^{4-d}}{i\pi^2} \int d^dk \frac{1; k_\mu; k_\mu k_\nu}{(k^2-m_1^2)((k+p)^2-m_2^2)}. \end{align} Considering the covariance under the Lorentz transformation, we can define another set of functions, \begin{align} B_\mu &\equiv p_\mu B_1,\\ B_{\mu\nu} &\equiv p_\mu p_\nu B_{21} + g_{\mu\nu} B_{22}. \end{align} There are three relations among them: \begin{align} p^2 B_1 &= \frac{1}{2}\left[ A(m_1) - A(m_2) + (-p^2 + m_2^2 - m_1^2) B_0 \right],\\ p^2 B_{21} + B_{22} &= \frac{1}{2} \left[ A(m_2) + (-p^2 + m_2^2 - m_1^2) B_1 \right],\\ p^2 B_{21} + d B_{22} &= A(m_2) + m_1^2 B_0, \end{align} through with all vector and tensor integrals are decomposed into linear combinations of scalar integrals $A$ and $B_0$. Thus, we only need to calculate the expression of the scalar integral $B_0$, which is given by \begin{align} B_0 &= \frac{1}{\epsilon} + 2 + x_1 \log \left( \frac{x_1-1}{x_1} \right) + x_2 \log \left( \frac{x_2-1}{x_2} \right) - \log (1-x_1) - \log (x_2-1) + \log \left( \frac{\tilde{\mu}^2}{-p^2} \right), \end{align} where $x_i$ are the roots of $-p^2 x^2 + (p^2 - m_1^2 + m_2^2) x + m_1^2 = 0$. It is also useful to see the massless cases \begin{align} B_0(p,m^2,0) &= \frac{1}{\epsilon} + 2 + \ln \frac{\tilde{\mu}^2}{-p^2+m^2} + \frac{m^2}{p^2}\ln\frac{-p^2+m^2}{m^2},\\ B_0(p,0,0) &= \frac{1}{\epsilon} + 2 + \ln \frac{\tilde{\mu}^2}{-p^2},\\ \end{align} and the large mass limit \begin{align} B_0(0,m_i^2,m_j^2) &= \frac{1}{\epsilon} + 1 + \frac{m_i^2}{m_i^2 - m_j^2} \ln \frac{\tilde{\mu}^2}{m_i^2} - \frac{m_j^2}{m_i^2 - m_j^2} \ln \frac{\tilde{\mu}^2}{m_j^2},\\ B_0(0,M^2,M^2) &= \frac{1}{\epsilon} + \ln \frac{\tilde{\mu}^2}{M^2},\\ B_0(0,0,M^2) &= \frac{1}{\epsilon} + 1 + \ln\frac{\tilde{\mu}^2}{M^2}, \end{align} where we assume $m_{i,j}$ with different indices take different non-zero values. Note that the PaVe functions in the large mass limit obey the recurrence relations \begin{align} B_0(0,m_i^2,m_j^2) &= \frac{1}{m_i^2 - m_j^2} \left[ A_0(m_i) - A_0(m_j) \right], \end{align} and similar for higher point functions.3-point functions
We define \begin{align} C_0(p_1^2,p_{12},p_2^2,m_1^2,m_2^2,m_3^2) \equiv \frac{\mu^{4-d}}{i\pi^2} \int d^dk \frac{1}{(k^2-m_1^2)((k+p_1)^2-m_2^2)((k+p_2)^2-m_3^2)}, \end{align} with $p_{12}\equiv p_1\cdot p_2$, and similar for vector and tensor integrals. The most general formula is too complicated to show and use, so we only focus on the specific limits here. By taking the large mass limit, we obtain \begin{align} C_0(0,0,0,m_i^2,m_j^2,m_k^2)&= \frac{m_i^2}{(m_i^2-m_j^2)(m_i^2-m_k^2)} \ln\frac{\tilde{\mu}^2}{m_i^2} + (\text{cycl.}),\\ C_0(0,0,0,m_i^2,m_i^2,m_j^2)&= -\frac{1}{m_i^2-m_j^2} + \frac{m_j^2}{(m_i^2-m_j^2)^2} \ln\frac{m_i^2}{m_j^2},\\ C_0(0,0,0,M^2,M^2,M^2)&=-\frac{1}{2M^2},\\ C_0(0,0,0,0,0,M^2)&=\frac{1}{M^2}\left[ \frac{1}{\epsilon_{\mathrm{IR}}} + \ln\frac{\tilde{\mu}^2}{M^2} + 1\right], \end{align} where $\epsilon_{\mathrm{IR}}=\frac{4-d}{2}$ is basically the same as $\epsilon$, but we add the subscript to clarify that it expresses the IR divergence in the limit $\epsilon_{\mathrm{IR}}\to 0$. We also obtain \begin{align} C_\mu(0,0,0,M^2,0,0)=-\frac{p_{1\mu} + p_{2\mu}}{M^2}\left[ \frac{1}{\epsilon_{\mathrm{IR}}} + \ln\frac{\tilde{\mu}^2}{M^2}\right]. \end{align}4-point functions
We define \begin{align} D_0(p_1^2,p_2^2,p_3^2,p_{12},p_{23},p_{31},m_1^2,m_2^2,m_3^2,m_4^2) \equiv \frac{\mu^{4-d}}{i\pi^2} \int d^dk \frac{1}{(k^2-m_1^2)((k+p_1)^2-m_2^2)((k+p_2)^2-m_3^2)((k+p_3)^2-m_4^2)}. \end{align} In the large mass limit, we obtain \begin{align} D_0(0,0,0,0,0,0,m_i^2,m_j^2,m_k^2,m_\ell^2)&= \frac{m_i^2}{(m_i^2-m_j^2)(m_i^2-m_k^2)(m_i^2-m_\ell^2)} \ln\frac{\tilde{\mu}^2}{m_i^2} + (\text{cycl.}),\\ D_0(0,0,0,0,0,0,m_i^2,m_i^2,m_j^2,m_j^2)&= -\frac{2}{(m_i^2-m_j^2)^2} + \frac{m_i^2+m_j^2}{(m_i^2-m_j^2)^3} \ln\frac{m_i^2}{m_j^2}.\\ \end{align}1. G. Passarino and M. Veltman, Nucl. Phys. B 160 (1979), 151-207 ↩