Hund's rule had been established empirically for the ground state of atoms and for configuration with equivalent electrons and it basically predicts the state term symbol $^{2S+1}L$ (Russel-Saunders notation). These terms must be introduced since every atomic state has many components that differ for the values of the orbital angular momentum $L$ and the intrinsic angular momentum $S$.
You may know that within the framework of the central field approximation the N-electron wave function is described by a Slater determinant and in general is a function of the spatial and spin coordinates. This antisymmetric wave function couples the spins of the electrons to the energy of the system.
Now I'm going to show a simplified example (I also made a list of warnings at the end) with 2 electrons in which there's an evidence of this coupling, this will help you get an idea of what the Hund's rule is based on. Consider a two electron atom described by an hamiltonian
$$H=-\frac{1}{2}\nabla^2_{\vec{r}_1}-\frac{1}{2}\nabla^2_{\vec{r}_2}-\frac{Z}{|\vec{r}_1|}-\frac{Z}{|\vec{r}_2|}+\frac{1}{|\vec{r}_{12}|}$$
and apply the perturbation theory in order to study the discrete exited states. We can treat the term $1/\vec{r}_{12}$ as a perturbation since it gives a small contribution to the energy and we can also see that the hamiltonian is spin independent, so the zero-order wave functions must only depend on the spatial coordinates $\Psi^{(0)}(\vec{r}_1, \vec{r}_2)$. We also have to take into account the exchange degeracy of the electron labels.
$$\Psi^{(0)}(\vec{r}_1, \vec{r}_2)=\frac{1}{\sqrt{2}}\left (\psi_{100}(\vec{r}_1)\psi_{nlm}(\vec{r}_2) \pm \psi_{nlm}(\vec{r}_1)\psi_{100}(\vec{r}_2)\right)$$
If we proceed with the calculation of the first perturbative term in energy using these two wave functions we obtain $E^{(1)}=J_{nl}\pm K_{nl}$, where the quantity $J_{nl}$ is known as direct or Coulomb integral and $K_{nl}$ is known as exchange integral, these terms are of the same order of magnitude. Using Pauli's matrices labelled with the numbers $1$ and $2$ standing for the electrons, one can observe that since we have $\sigma_1 \cdot \sigma_2=-3\mathbb{1}_{d}$ for a spin singlet and $\sigma_1 \cdot \sigma_2=\mathbb{1}_{d}$ for a triplet spin state, we can rewrite the energy correction as follows:
$$E^{(1)}=J_{nl}-\frac{1}{2}(\mathbb{1}_{d}+\sigma_1 \cdot \sigma_2) K_{nl}$$
As a consequence, singlet states always have a positive correction while the triplet states always bring a negative correction, this conclusion can also be generalized to more electron atoms. The Pauli exclusion principle introduces a coupling between the space and spin variables and the total effect is that electrons move under the influence of a force whose sign depends on the relative orientation of the spin. Remember that the singlet state correspond to zero unpaired electrons and the triplet correspond to two unpaired electrons. Now for example if you have 2 electrons and the orbitals $2p_{x}, 2p_{y}, 2p_{z}$ to fill, if you think in terms of stability of the atom the electrons will occupy the state with lower energy, so they will be unpaired.
Warnings:
The trend of an atom in absence of EM field is to stay at the lowest possible energy level, so we will normally find the system at the ground state. An atom with 2 electrons will therefore stay in the $1S^{2}$ orbital. Despite this fact, I decided to show the case of the energy corrections to the excited levels since it is an example of the typical coupling spin-space which occurs also in the other systems with more electrons.
In this particular case we have assumed that $1/|\vec{r}_{12}|$ was a perturbation but for atoms with more than 2 electrons we can't do this approximation. To solve the problem for atoms with more electrons we have to deal with an effective potential which takes into account the attraction of the nucleus and the repulsive effect of the electrons, the Hartree-Fock potential satisfies these conditions.
