If we solve the time-independent Schrödinger equation for any atom by considering only the electrostatic potential, an electron has the same probability of occupying the $p_{x}$ orbital as it does the $p_{y}$ and $p_{z}$ orbitals, because the electrostatic force is a central / spherically symmetric force which depends only on $r$, and the wavefunctions of the $p_{x}$ , $p_{y}$ , $p_{z}$ orbitals only differ in terms of the angular coordinates $(\theta,\phi)$.
However, a nucleus has a quantum spin which introduces a magnetic field, and the symmetry is broken. If we take into account the magnetic field of the nucleus, then the electromagnetic force which the nucleus exerts on electrons occupying the $p_{x}$, $p_{y}$, and $p_{z}$ orbitals is no longer central. So there should be a difference between the energies of the wavefunctions of the $p_{x}$, $p_{y}$, and $p_{z}$ orbitals.
How significant is this effect? In atoms with a small atomic number or with a small nucleus spin the effect is negligible and can be only observed near absolute zero. However shouldn't it become significant for atoms with bigger atomic number or with a large nuclear spin?