Any suggestions how to get started? I know Liouville's theorem, but not sure how to apply it here:
Let $u$ be a harmonic function on $\mathbb{R}^3$. Assume there exists $C>0$, independent of $x$, such that $$|u(x)| \leq C(1+|x|)$$ on $\mathbb{R}^3$. Show that then $\partial_i u$ is constant on $\mathbb{R}^3$, $\forall i = 1,2,3$
First, balls versus spheres: The ball with center at the origin and radius $r$ is $$B_r=\{x\,:\,|x|<r\},$$as opposed to the sphere of radius $r$, defined by $|x|=r$.
We know that harmonic functions satisfy the mean-value property in the sense that $u(x)$ is the average of $u$ over a sphere centered at $x$. If you integrate in "polar coordinates" you get a mean-value property for balls:
$$u(x)=\frac{c}{r^3}\int_{x+B_r}u.$$Here $x+A=\{x+y\,:\,y\in A\}$, so that $x+B_r$ is the ball of radius $r$ centered at $x$.
Given two sets $A$ and $B$, the symmetric difference is $$A\Delta B=(A\setminus B)\cup(B\setminus A).$$That is, $A\Delta B$ is the set of points in one of the sets but not both.
Now say $x\in\Bbb R^3$ and take $r>|x|$. To be on the safe side say $r>2(1+|x|)$. Now $$u(0)-u(x)=\int_{B_r}u-\int_{x+B_r}u$$. So $$|u(x)-u(0)|\le\frac c{r^3}\int_{B_r\Delta(x+B_r)}|u|\le\frac c{r^3}M_{r+|x|}m(B_r\Delta(x+B_r)),$$where $m(S)$ is the volume of $S$ and $M_r$ is the supremum of $|u(y)|$ for $|y|\le r$.
The hypothesis on $u$ shows that $M_{|x|+r}\le 1+|x|+r$. If $r$ is large enough this shows $$M_{|x|+r}\le 2r.$$
Now note that $x+B_r\subset B_{|x|+r}$. Hence $$(x+B_r)\setminus B_r\subset B_{|x|+r}\setminus B_r.$$ That's a spherical shell with radius $|x|+r$ and thickness $|x|$, so again if $r$ is much larger than $x$ we get $$m((x+B_r)\setminus B_r)\le cr^2|x|$$(note the standard convention in analysis: the value of "$c$" cah change from line to line...) Similarly for the other half of $B_r\Delta(x+B_r)$, so $$m((x+B_r)\Delta B_r)\le cr^2|x|.$$
Put all these inequalities together and divide by $|x|$ and you get $$\frac{|u(0)-u(x)|}{|x|}\le c.$$This shows that $$|\partial_ju(0)|\le c.$$
The same argument gives the same bound at every other point. So the partials of $u$ are bounded, hence constant.
(If you work out the details it may look like you get a different bound at other points. Let $r\to\infty$.)