First, we solve it for the unit sphere, since the solution is just scaled up by $a$.
Secondly, we observe that if we have a single octant, with center of mass at $(u, u, u)$, then if we combine the four positive-$z$ octants (say), then the center of mass will be at $(0, 0, u)$, by symmetry. So if we find the center of mass of the positive-$z$ hemisphere, we thereby obtain the center of mass of any octant.
We note that the center of mass of a unit hemispherical shell is at coordinate
$$
\frac{\int_{\theta=0}^{\pi/2} \cos\theta\sin\theta \, d\theta}
{\int_{\theta=0}^{\pi/2} \cos\theta \,d\theta} = \frac{1}{2}
$$
so that of a hemispherical shell of radius $r$ is at coordinate $\frac{r}{2}$.
In the sphere in question, we have shells of area proportional to $r^2$, multiplied by a density proportional to $r^3$, so the coordinate of the center of mass of the hemisphere is
$$
\frac{\int_{r=0}^1 \frac{r}{2} r^5 \, dr}{\int_{r=0}^1 r^5 \, dr}
= \frac{\frac{1}{14}}{\frac{1}{6}} = \frac{3}{7}
$$
For a uniform sphere, it would be
$$
\frac{\int_{r=0}^1 \frac{r}{2} r^2 \, dr}{\int_{r=0}^1 r^2 \, dr}
= \frac{\frac{1}{8}}{\frac{1}{3}} = \frac{3}{8}
$$
thus confirming the values obtained by Ron Gordon. It is mildly curious to me that that high a density gradient has such a small effect.
$$\bar{\cos{\theta}} = \int dr \, r^2 f(r) \, \int d\phi \, \sin{\phi} \, \int d\theta \, \cos{\theta}$$
$$\bar{\sin{\phi}} = \int dr \, r^3 f(r) \, \int d\phi \, \sin^2{\phi} \, \int d\theta \, $$
$$\bar{r} = \int dr \, r^3 f(r) \, \int d\phi \, \sin{\phi} \, \int d\theta \, $$
Then
$$\begin{align}\bar{r} \,\bar{\sin{\phi}}\,\bar{\cos{\theta}} &= \frac{\int dr \, r^3 f(r)\left (\int dr \, r^3 f(r) \right )^2 \int d\phi \, \sin^2{\phi} \left ( \int d\phi \, \sin{\phi} \right )^2\int d\theta \, \cos{\theta} \left ( \int d\theta \right )^2}{\left (\int dr \, r^2 f(r) \right )^3 \left ( \int d\phi \, \sin{\phi} \right )^3 \left ( \int d\theta \right )^3}\\ &= \frac{\int dr \, r^3 f(r) \, \int d\phi \, \sin^2{\phi} \,\, \int d\theta \, \cos{\theta} }{\int dr \, r^2 f(r)\, \int d\phi \, \sin{\phi}\, \int d\theta }\\ &= \bar{x}\end{align}$$
But I do not think you can work with the bare $\phi$ and $\theta$ because the trig functions introduce a nonlinearity that belies the linearity of the integrals.
Best Answer
Center of mass calculation is essentially the same as other averages such as expectation value in probability theory. Essentially it is the weighted average of the space coordinate.
The definition is best wrote down in Cartesian coordinates:
$$\mathbf{r}_\text{center} = \int_V \mathbf{r} f(\mathbf{r}) \mathrm{d} V$$
with $\mathbf{r} = (x,y,z)$ and $\mathrm{d} V = \mathrm{d} x \mathrm{d} y \mathrm{d} z$. And $$f(\mathbf{r}) = \frac{1}{\int_V\rho(\mathbf{r}) \mathrm{d} V} \rho(\mathbf{r})$$ the volume normalized distribution function.
This hint should allow you to solve the stated problem. It is now just a matter of correct calculation of the integrals with spherical coordinates. You will get the center of mass in Cartesian coordinates, though. But I dont if there exists any simple integral expression to obtain the center of mass directly in spherical coordinates.