What do moments of inertia do in the potential terms of Lagrangians? I am struggling to understand the Lagrangian computed in this paper. In particular, a binary spacecraft-debris system is assumed as below.

The analysis goes as follows.

1- I am in trouble to understand the terms highlighted in the second image above. Namely, I had never seen any term including moments of inertia in potential elements of any Lagrangian. What are the meanings of those terms? In particular, I can somehow grasp the potential nature of the red term. Given
$$I_z = I_x + I_y$$
That term becomes
$$- \dfrac{\mu(I_x + I_y+I_z)}{2r_{B}^{3}} = - \dfrac{\mu(2m_{B}r_{B}^{2})}{2r_{B}^{3}} = - \dfrac{\mu m_{B}}{r_{B}},$$
Which is the same as the second term of equation (3)! Can you please explain why would one consider this term and what it is called? 
2- The blue term is totally cryptic to me! I cannot resolve not only the meaning of this term but also its sign (which is positive contrary to other potential terms). Any comments on how it has been derived are appreciated.
 A: The gravitational potential energy of a material point is given by:
$$
U=-\mu\frac{m_0}{R}.
$$
However, most of the bodies are not material points. Thus, the potential energy can be given by an integral:
$$
U=-\mu\int\frac{dm}{r}=-\mu\int\frac{dm}{|\pmb R+\pmb x|},
$$
where $R$ is the position of the centre of mass (CoM) and $x$ is the coordinate of point on CoM coordinate system. Series expansion gives us:
$$
\frac{1}{|\pmb R+\pmb x|}=\frac1{R}-\frac{R_\alpha}{R^3}x_\alpha-\frac{3R_\alpha R_\beta-\delta_{\alpha\beta}R^2}{R^5}\frac{x_\alpha x_\beta}2 +\ldots
$$
(Einstein summation assumed)
The integral of the first term gives the potential of a material point. The integral of the second term vanishes, because $R$ is CoM and $\int x_\alpha dm=0$. The integral of the third one:
$$
-\frac{\mu}{2R^5}\int(3R_\alpha R_\beta-\delta_{\alpha\beta}R^2)x_\alpha x_\beta dm =\\ -\frac{\mu R_\alpha R_\beta}{2R^5}\int(3x_\alpha x_\beta-\delta_{\alpha\beta}x^2)dm = -\frac{\mu}{2R^5}Q_{\alpha\beta}R_\alpha R_\beta
$$
Tensor $Q_{\alpha\beta}$ is called gravitational quadrupole moment, and the technique we used is called multinomial expansion.
Recall that the tensor of inertia:
$$
I_{\alpha\beta}=\int(\delta_{\alpha\beta}x^2-x_\alpha x_\beta)dm.
$$
Then one can show that $Q_{\alpha\beta} = I_{\gamma\gamma}\delta_{\alpha\beta}-3I_{\alpha\beta}$.
Thus the quadrupole term of potential energy:
$$
U^{(2)} = -\mu\frac{\mathop{\mathrm{Tr}}I}{2R^3}+3\mu\frac{I_{\alpha\beta}n_\alpha n_\beta}{2R^3},
$$
where $n_\alpha=R_\alpha/R$ is a unit vector. The first term is exactly what you see in the paper, the second one can be likely transformed to the desired form with the help of trigonometry.
By the way, the problem with your formula in (1) is that $r_B$ is the position of the body relative to the earth, whereas moment of inertia $I$ is given about CoM, so $\mathop{\mathrm{Tr}}I\neq 2m_Br_B^2$.
