Where do symmetries in atomic orbitals come from?

The hydrogen atom is spherically symmetric, so for any solution of the Schrödinger equation for the hydrogen atom, any rotation of that solution must also be a solution. If you do the math on how to rotate a solution, it turns out that the solutions with a particular energy $E_n$ fall into groups labeled by an integer $l < n$. The integer $l$ is physical: $\hbar^2 l(l+1)$ is the magnitude squared of the angular momentum. Within each group, rotating the solution gives you a new solution in the same group. These two facts are of course connected: a rotation can't change the length of a vector.

One can show that each group contains $2l+1$ independent solutions, in that any solution $|n,l\rangle$ where the energy is $E_n$ and the angular momentum $\hbar^2 l(l+1)$ can be written as a sum $$|n,l\rangle = \sum_{m=-l}^l c_m |n,l,m\rangle$$ (I apologize for the somewhat poor notation.)

This decomposition is based on choosing a particular axis, and taking each state to depend on the angle $\varphi$ around this axis as $e^{im\varphi}$. The appearance of axes of symmetry in these plots is due to this choice of axis and particular decomposition. With another choice of axis, which is the same as a rotation, the states will be mixed.

The bottom line is that it's not each solution -- wavefunction -- that needs to be spherically symmetric, but the total set of solutions.

How do these symmetries shown in the above article occur? What about the 'preferable' axis of symmetries? Why these?

For atoms subject to no net external electric of magnetic fields the orientation of the axes is arbitrary. This shows up clearly in the math because adding up all the spherical harmonic contributing to a single shell (1s, 2s, 2p, 3s, 3p, 3d, ...) gives no angular dependence. It doesn't show clearly in the visualization because those plots employ an arbitrary cut-off in generating the display. So, short answer, the lobes of the orbitals point along the coordinate axes purely for convenience: there is no physics content to that feature of the rendering.

The fact that there are a non-negative integer number of radial or angular nodes arises from the boundary conditions on the wave-function: just like the vibrations of a guitar string only those modes that 'fit' in the space exist as time-independent solutions.

In the case that there are external electromagnetic fields, then those fields do two things:

• They change the shape of the time-independent solutions
• The enforce a choice of orientation on the new solutions.