If we have the Laplace equation for electrostatics in free space, that is
$$\Delta u(x) = 0 \quad \quad x \in \mathbb{R}^3,$$
is the only solution $u = 0$? And also, we only get non-zero solutions for $u$ if we instead consider the Laplace equation on some bounded domain? How can I mathematically show this if this is in fact the case?
Since you are dealing with the whole space, you can take advantage of the so called Liouville theorem for harmonic functions.
Liouville theorem. Let $\phi : \mathbb R^n \to \mathbb R$ be a $C^2$ function such that $\Delta \phi =0$ everywhere. If $\phi$ is bounded (i.e., there is $k \in [0,+\infty)$ such that $|\phi(x)| \leq k$ for every $x \in \mathbb R^n$), then $\phi$ is everywhere constant.
This has a pair of corollaries.
Corollary 1. Let $\phi : \mathbb R^n \to \mathbb R$ be a $C^2$ function such that $\Delta \phi =0$ everywhere. If $\phi$ is bounded and $\phi(a n) \to 0$ as $a \to +\infty$, where $n \in \mathbb R^n$ is a fixed unit vector, then $\phi=0$.
Proof. $\phi$ is bounded, thus $\phi(x)=c$ constantly due to Liouville theorem. $c = \lim_{a\to +\infty} c = \lim_{a\to +\infty} \phi(a n) = 0$.
Corollary 2. Let $\phi : \mathbb R^n \to \mathbb R$ be a $C^2$ function such that $\Delta \phi =0$ everywhere. If $\phi(x)$ tends to $0$ as $|x|\to +\infty$ (i.e., for every $\epsilon >0$ there is $r_\epsilon>0$ such that $|\phi(x)|< \epsilon$ if $|x|> r_\epsilon$), then $\phi=0$.
Proof. Take $\epsilon >0$ so that $|\phi(x)|< \epsilon$ if $|x|> r_\epsilon$. In the compact set $|x| \leq 2r_\epsilon$, $\phi$ is continuous (since it is $C^2$) and thus it is bounded therein by some $M \geq 0$. Consequently $|\phi(x)| \leq \epsilon_r + M$ for all $x \in \mathbb R^n$. Liouville theorem now implies that $\phi(x)=c$ constantly. However this constant $c$ must satisfy $0 \leq |c |< \epsilon$ for every $\epsilon >0$ and thus $c=0$.
The second corollary uses a very weak requirement regarding how $\phi$ uniformly tends to $0$ for $|x|\to +\infty$. Obviously $\phi(x) \sim const/|x|$ is OK, but also much weaker convergences are enough, like $ \sim const / \ln |x|$...
I'm not sure that I understand what you mean by "bounded domain". This is an answer based on what I understand from the question.
To solve Laplace's equation uniquely, you have to specify either the Dirichlet or Neumann type boundary condition.
EDIT: The boundary of the region of interest can be at finite distances as well as at infinity. For example, to solve the potential due to a point change in free space, you have to solve Poisson's equation. Here, we take the potential going to zero at infinity. So if you have some finite boundary in mind that is not necessary. As you know the solution in this case is the familiar $u(r)\sim \frac{1}{r}$ Coulomb potential (which is not identically zero everywhere).
Now for Laplace's equation in absolutely free space (no charge anywhere), if the boundary condition is such that the potential vanishes everywhere on the boundary, then the potential will remain zero everywhere simply because Laplace's equation doesn't support local maxima or minima.
