This is a community-wiki answer trying to remove this question from the unanswered queue.
I am not that comfortable with indices notation, so I will translate the integral into traditional multivariate calculus language. Your $\vec{x}$ will be bold face letter $\mathbf{r} = (x,y,z)$, $\mathbf{a} = (a_x,a_y,a_z)$, and $x_i x_j$ will be $xy$ in my answer if $i=1$, and $j=2$. The subscript $(\cdot)_z$ of a vector will just be its $z$-component.
$2\int d^3x \,x_i a_j(\vec x)=\epsilon_{ijk} \Big[\int \,d^3x \,\, \vec x\times \vec a(\vec x)\Big]_k -i b \int\,d^3x\,\, x_i x_j c(\vec x)$
Knowing $-ibc = \nabla \cdot \mathbf{a}$. This integral can be written as six integral identities if counting all the permutation and get rid of the Levi-Civita symbol: $ 2\int_{\mathbb{R}^3} x a_y \,dxdydz = \int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_z \,dxdydz - \int_{\mathbb{R}^3} xy\nabla \cdot \mathbf{a}\,dxdydz, \\ 2\int_{\mathbb{R}^3} y a_x \,dxdydz = -\int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_z \,dxdydz - \int_{\mathbb{R}^3} yx\nabla \cdot \mathbf{a}\,dxdydz, \\ 2\int_{\mathbb{R}^3} y a_z \,dxdydz = \int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_x \,dxdydz - \int_{\mathbb{R}^3} yz\nabla \cdot \mathbf{a}\,dxdydz, \\ 2\int_{\mathbb{R}^3} z a_y \,dxdydz = -\int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_x \,dxdydz - \int_{\mathbb{R}^3} zy\nabla \cdot \mathbf{a}\,dxdydz, \\ 2\int_{\mathbb{R}^3} z a_x \,dxdydz = \int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_x \,dxdydz - \int_{\mathbb{R}^3} zx\nabla \cdot \mathbf{a}\,dxdydz, \\ 2\int_{\mathbb{R}^3} x a_z \,dxdydz = -\int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_x \,dxdydz - \int_{\mathbb{R}^3} xz\nabla \cdot \mathbf{a}\,dxdydz. $ Here I will check if the first one is correct. Others can be checked in the exactly same fashion. Now $\mathbf{r}\times \mathbf{a} = (ya_z - za_y, za_x - xa_z, xa_y-ya_x).$ The right hand side of the first identity is: $ \text{R.H.S}= \int_{\mathbb{R}^3} (\mathbf{r}\times \mathbf{a})_z \,dxdydz - \int_{\mathbb{R}^3} xy\nabla \cdot \mathbf{a}\,dxdydz \\ = \int_{\mathbb{R}^3} (xa_y-ya_x) \,dxdydz + \int_{\mathbb{R}^3} \nabla(xy)\cdot \mathbf{a}\,dxdydz, $ in which we assume the decaying condition of $|\mathbf{a}| = o(r^{-4})$ where $r = |\mathbf{r}|$, this assumption just guarantees that if we are integrating on a really large ball in $\mathbb{R}^3$, the boundary integral term produced by Gauss-Divergence theorem vanishes when the radius of the ball goes to infinity. Expanding we have: $ \text{R.H.S} = \int_{\mathbb{R}^3} (xa_y-ya_x) \,dxdydz + \int_{\mathbb{R}^3} (y,x,0)\cdot \mathbf{a}\,dxdydz \\ =\int_{\mathbb{R}^3} (xa_y-ya_x) \,dxdydz + \int_{\mathbb{R}^3} (ya_x+xa_y)\,dxdydz \\ = 2\int_{\mathbb{R}^3} x a_y \,dxdydz = \text{L.H.S}. $