Chemistry - Why is there a z² label for d orbitals, but no x² and y² labels with corresponding shapes?

TL;DR The $\mathrm{d}_{z^2}$ orbital is a result of solving the Schrödinger equation for the hydrogen atom in the most mathematically convenient way.

To properly understand this, it is necessary to go back to the fundamentals. The complex d-orbitals are obtained by solution of the Schrödinger equation. In general, these d-orbitals are made up of a radial part $R(r)$ and an angular part $Y(\theta,\phi)$. The radial part is not relevant to the question, because its value does not depend on the direction in which you look, so it cannot explain why there is a $\mathrm{d}_{z^2}$, but not a $\mathrm{d}_{x^2}$. So, it is instructive to look at the angular component (the spherical harmonics), which are tabulated below. The normalisation factors are ignored here:

$$\begin{align} \mathrm{d}_{+2} &\propto R(r) \cdot \sin^2\theta \cdot \mathrm{e}^{\mathrm 2i \phi} \\ \mathrm{d}_{+1} &\propto R(r) \cdot (-\sin\theta \cos\theta) \cdot \mathrm{e}^{\mathrm i \phi} \\ \mathrm{d}_0 &\propto R(r) \cdot (3\cos^2\theta - 1) \\ \mathrm{d}_{-1} &\propto R(r) \cdot \sin\theta \cos\theta \cdot \mathrm{e}^{-\mathrm i \phi} \\ \mathrm{d}_{-2} &\propto R(r) \cdot \sin^2\theta \cdot \mathrm{e}^{\mathrm -2i \phi} \end{align}$$

These are combined in an appropriate manner to generate the real d-orbitals. It turns out that the $\mathrm{d}_{z^2}$ orbital is simply the same as the original $\mathrm{d}_0$ orbital, which is already real (the radial component $R(r)$ is real). Since $\cos\theta = z/r$ and $r^2 = x^2 + y^2 + z^2$ (spherical coordinates), the angular component can be rewritten as

$$\begin{align} 3\cos^2\theta - 1 &= 3\left(\frac{z}{r}\right)^2 - 1 \\ &= \frac{3z^2 - r^2}{r^2} \\ &= \frac{2z^2 - x^2 - y^2}{r^2} \end{align}$$

and conventionally we denote this as $\mathrm{d}_{z^2}$. There is no analogous $\mathrm{d}_{x^2}$ orbital because the conventional way of deriving real d-orbitals does not yield any orbital with an angular form of $(2x^2 - y^2 - z^2)/r^2$.

But that still does not properly answer the question. To do so, we should look more closely at how the complex d-orbitals are obtained. The Schrödinger equation, $\hat{H}\psi = E\psi$, only prescribes that we need to find eigenstates of $\hat{H}$, but does not say how. It turns out that the easiest way to accomplish this is to look for solutions which are simultaneous eigenstates of $\hat{H}$ (energy), $\hat{L}^2$ (total orbital angular momentum), and $\hat{L}_z$ (projection of orbital angular momentum onto the z-axis). The eigenvalues of the solutions under these operators are directly tied to the quantum numbers $n$, $l$, and $m_l$, respectively.

Fundamentally, the $\mathrm{d}_{0}$ orbital (i.e. the $\mathrm{d}_{z^2}$ orbital) arises because we have chosen it to be an eigenstate of the operator $\hat{L}_z$. If we were to solve the Schrödinger equation in a different way, by searching for simultaneous eigenstates of $\hat{H}$, $\hat{L}^2$, and $\hat{L}_x$ (instead of $\hat{L}_z$), we would indeed get a $\mathrm{d}_{x^2}$ orbital. This is an equally valid solution of the Schrödinger equation (as it is an eigenstate of $\hat{H}$), but it is simply not the conventional solution.

So, why choose $\hat{L}_z$ and not $\hat{L}_x$ when solving the Schrödinger equation? The answer is simply mathematical convenience: the Schrödinger equation is most easily solved in spherical coordinates, and in spherical coordinates, $\hat{L}_z$ has a particularly nice form of being equal to $-\mathrm{i}\hbar(\partial/\partial \phi)$. The eigenstates of this are none other than $\mathrm{e}^{\mathrm{i}m\phi}$, which is exactly what we observe in the complex d-orbitals at the beginning of this post ($m$ is, of course, quantised, so there are only five valid solutions with $l = 2$, i.e. d-orbitals).

But why does $\hat{L}_z$ have a nice form and not $\hat{L}_x$? Ultimately, this mathematical convenience is a result of the definition of spherical coordinates. Because $\theta$ is chosen to be the angle made by the vector with the z-axis, the z-axis has a special status in spherical coordinates:

$$\begin{align} z &= r\cos\theta \\ x &= r\sin\theta\cos\phi \\ y &= r\sin\theta\sin\phi \end{align}$$

So, at the end of the day, it is a completely arbitrary convention. And it should be, too, because no axis should have an inherently special status in a spherically symmetric system like the hydrogen atom. As mentioned previously, you could obtain a $\mathrm{d}_{x^2}$ orbital by solving the Schrödinger equation in a different way. You could even choose to redefine the axes in such a way which gave you a $\mathrm{d}_{x^2}$ orbital (it's quite easy; simply rename $(x,y,z)$ to $(y,z,x)$).

But there is no way you can have both at the same time. Because $\hat{L}_x$ and $\hat{L}_z$ do not commute with one another, it is not possible to "solve" the Schrödinger equation in a way which gives you both $\mathrm{d}_{x^2}$ and $\mathrm{d}_{z^2}$ orbitals at the same time.