What is the physical basis for Hund's first rule?
Disclaimer

My following answer is the "traditional" explanation of Hund's first rule, which is based on a smaller value of $V_\mathrm{ee}$ (electron-electron repulsions) in the triplet state arising from Fermi holes. According to Levine's Quantum Chemistry 7th ed.:


  This traditional explanation turns out to be wrong in most cases. It is true that the probability that the two electrons are very close together is smaller for the helium $^3S$ $1s2s$ term than for the $^1S$ $1s2s$ term. However, calculations with accurate wave functions show that the probability that the two electrons are very far apart is also less for the $^3S$ term. The net result is that the average distance between the two electrons is slightly less for the $^3S$ term than for the $^1S$ term, and the interelectronic repulsion is slightly greater for the $^3S$ term.
  
  The calculations show that the $^3S$ term lies below the $^1S$ term because of a substantially greater electron nucleus-attraction in the $^3S$ term as compared with the $^1S$ term. Similar results are found for terms of the atoms beryllium and carbon.


Wildcat has an answer which briefly mentions this. Aside from that Levine gives several references:


Katriel, J.; Pauncz, R. Theoretical Interpretation of Hund's Rule. Adv. Quantum Chem. 1977, 10, 143–185. DOI: 10.1016/S0065-3276(08)60580-8.
Shim, I.; Dahl, J. P. Theor. Chim. Acta 1978, 48 (2), 165. The DOI doesn't resolve properly, but here's the link to the publisher's website.
Boyd, R. J. A quantum mechanical explanation for Hund's multiplicity rule. Nature 1984, 310, 480–481. DOI: 10.1038/310480a0.
Oyamada, T.; Hongo, K.; Kawazoe, Y.; Yasuhara, H. Unified interpretation of Hund's first and second rules for 2p and 3p atoms. J. Chem. Phys. 2010, 133, 164113. DOI: 10.1063/1.3488099.


With that said, here is the old-fashioned argument for why parallel spins are favoured over paired spins (all else being equal).



Introduction

The origin of Hund's first rule lies in "exchange energy", which is a way of saying that electrons with like spin repel each other less than electrons with unlike spin.

To see why, let's consider an excited state of the helium atom: $\mathrm{1s^1 2s^1}$. (The ground state configuration $\mathrm{1s^2}$ is of no use to us because those electrons have to be paired with unlike spin.) In the excited state, the electrons can either have parallel spin (the triplet case) or paired spin (the singlet case).

The total electronic wavefunction comprises a spatial part, which describes which orbitals the electrons are in (i.e. it specifies $n,l,m_l$), and a spin part, which describes the, well, spin ($m_s$). According to the Pauli exclusion principle, the total wavefunction must be antisymmetric with respect to interchange of the electron labels: $\Psi(1,2) = -\Psi(2,1)$.

Spin wavefunctions

The spin wavefunctions are simultaneous eigenfunctions of the $\hat{S}^2$ and $\hat{S}_z$ operators. The total spin quantum number, $S$, can be found from the eigenvalue of the operator $\hat{S}^2$:

$$\hat{S}^2\psi = \hbar^2 S(S+1)\psi$$

For more detailed discussion, refer to any QM text on the coupling of angular momenta.

There are three possible triplet spin wavefunctions, which are degenerate.

$$\psi_{\mathrm{spin,triplet}}(1,2) = \begin{cases}
\alpha(1)\alpha(2) \\
\beta(1)\beta(2) \\
\frac{1}{\sqrt{2}}[\alpha(1)\beta(2) + \alpha(2)\beta(1)]
\end{cases}$$

It turns out that $\hat{S}^2\psi_\mathrm{triplet} = 2\hbar^2\psi_{\mathrm{triplet}}$, which means that these have the quantum number $S = 1$. This corresponds to a multiplicity of $M = 2S + 1 = 3$, hence "triplet". These wavefunctions represent the case where the electrons have parallel spin. (Loosely speaking, each electron has a spin of $1/2$, and if these are aligned parallel the total spin is $1/2 + 1/2 = 1$.)

The singlet state is

$$\psi_{\mathrm{spin,singlet}}(1,2) = \frac{1}{\sqrt{2}}[\alpha(1)\beta(2) - \alpha(2)\beta(1)]$$

and as you might guess, $\hat{S}^2\psi_\mathrm{singlet} = 0$, meaning that $S = 0$ and $M = 2S+1 = 1$. This represents the case where the electrons are paired - one can think of it as the two $1/2$ electrons spins “cancelling each other out”.

The important thing to notice is that the triplet spin wavefunctions are symmetric with respect to permutation of the labels, whereas the singlet spin wavefunction is antisymmetric. For example:

$$\begin{align}
\psi_\mathrm{singlet}(2,1) &amp;= \frac{1}{\sqrt{2}}[\alpha(2)\beta(1) - \alpha(1)\beta(2)] \\
&amp;= -\frac{1}{\sqrt{2}}[\alpha(1)\beta(2) - \alpha(2)\beta(1)] \\
&amp;= -\psi_\mathrm{singlet}(1,2)
\end{align}$$

Spatial wavefunctions

The only way of constructing these are the linear combinations below.

$$\begin{align}
\psi_\mathrm{space,symm} &amp;= \frac{1}{\sqrt{2}}[\mathrm{1s}(1)\mathrm{2s}(2) + \mathrm{1s}(2)\mathrm{2s}(1)] \\
\psi_\mathrm{space,antisymm} &amp;= \frac{1}{\sqrt{2}}[\mathrm{1s}(1)\mathrm{2s}(2) - \mathrm{1s}(2)\mathrm{2s}(1)]
\end{align}$$

Remember that the atomic orbitals above are functions of the coordinates of the electron: $\mathrm{1s}(1) = \mathrm{1s}(\vec{r}_1) = \mathrm{1s}(r_1,\theta_1, \phi_1)$. Now, the antisymmetric spatial wavefunction has an interesting property. Let's explore what happens when the electrons get very close together, i.e. $\vec{r}_1 = \vec{r}_2 = \vec{r}$.

$$\begin{align}
\psi_\mathrm{space,antisymm} &amp;= \frac{1}{\sqrt{2}}[\mathrm{1s}(\vec{r})\mathrm{2s}(\vec{r}) - \mathrm{1s}(\vec{r})\mathrm{2s}(\vec{r})] \\
&amp;= 0
\end{align}$$

If the wavefunction is 0 when the two electrons are on top of each other, that means that there is no probability of this happening. The average distance between the two electrons $\left\langle\left|\vec{r}_2 - \vec{r}_1\right|\right\rangle$ will therefore be larger, and on average the electrons will repel each other less, since they are on average further away from each other. This is known as a Fermi hole.

In fact, one can also show that the symmetric wavefunction has a larger probability amplitude when $\vec{r}_1 = \vec{r}_2$. This is known as a Fermi heap.

The total wavefunction

As I mentioned earlier, the Pauli exclusion principle dictates that the total wavefunction be antisymmetric. This means that the triplet spin wavefunctions, which are symmetric, must be paired with the antisymmetric spatial wavefunction. Likewise, the singlet spin wavefunction, which is antisymmetric, must be paired with the symmetric spatial wavefunction.

$$\begin{array}{c|c|c}
\text{Spin} &amp; \text{Spatial} &amp; \text{Total}\\
\hline
\text{Triplet (S)} &amp; \text{AS} &amp; \text{S } \times \text{ AS} = \text{AS} \\
\text{Singlet (AS)} &amp; \text{S} &amp; \text{AS } \times \text{ S} = \text{AS}
\end{array}$$

So, there is a decrease in electron-electron repulsion when the electrons are parallel - because this means they have a triplet spin wavefunction, and hence an antisymmetric spatial wavefunction, which has a Fermi hole.



Exchange integrals

There is another, perhaps more mathematical, way of looking at it, and that is to calculate the expectation value of the electronic energy for both wavefunctions. For the helium atom, under the assumption of an infinitely heavy nucleus (essentially an atomic version of the Born-Oppenheimer approximation), the total Hamiltonian (in atomic units) is

$$\hat{H} = -\frac{1}{2}(\nabla_1^2+\nabla_2^2) - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{r_{12}}$$

where $r_i$ is the distance of electron $i$ from the nucleus, and $r_{12}$ is the distance between the two electrons. The expectation value of the energy is simply $\left\langle\Psi\middle|\hat{H}\middle|\Psi\right\rangle$ (assuming our wavefunction is already normalised). It turns out that the only important term in our present discussion is the electron-electron repulsion term $1/r_{12}$.

Let's look at the triplet spin wavefunction (of course paired with the antisymmetric spatial wavefunction) first. The operator doesn't depend on spin degrees of freedom, so we can separate those, and the inner product of the spin wavefunction with itself is just unity.

$$\begin{align}
\left\langle\Psi_\mathrm{triplet}\middle|\frac{1}{r_{12}}\middle|\Psi_\mathrm{triplet}\right\rangle 
&amp;= \left\langle\psi_\mathrm{spin,triplet}\middle|\psi_\mathrm{spin,triplet}\right\rangle\left\langle\psi_\mathrm{space,antisymm}\middle|\frac{1}{r_{12}}\middle|\psi_\mathrm{space,antisymm}\right\rangle \\
&amp;= \frac{1}{2} \left\langle \mathrm{1s}(1)\mathrm{2s}(2) - \mathrm{1s}(2)\mathrm{2s}(1) \middle| \frac{1}{r_{12}} \middle| \mathrm{1s}(1)\mathrm{2s}(2) - \mathrm{1s}(2)\mathrm{2s}(1)  \right\rangle \\
&amp;= 
\frac{1}{2} \Biggr\{ \left\langle \mathrm{1s}(1)\mathrm{2s}(2) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(1)\mathrm{2s}(2) \right\rangle  \\
&amp;\quad\quad\quad - \left\langle \mathrm{1s}(2)\mathrm{2s}(1) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(1)\mathrm{2s}(2) \right\rangle \\
&amp;\quad\quad\quad - \left\langle \mathrm{1s}(1)\mathrm{2s}(2) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(2)\mathrm{2s}(1) \right\rangle \\
&amp;\quad\quad\quad + \left\langle \mathrm{1s}(2)\mathrm{2s}(1) \middle|\frac{1}{r_{12}}\middle| \mathrm{1s}(2)\mathrm{2s}(1) \right\rangle\Biggr\}
\end{align}$$

Ugly, but that's what it is.

We can simplify it however by noting that the first and the fourth brackets are the same (since the variables in the integral are just dummy variables, interchanging the electron labels does not affect the value of the integral). This integral is traditionally called a Coulomb integral, since it simply reflects the average Coulombic repulsion between two independent (uncorrelated) electrons, one in the 1s orbital and the second in the 2s orbital. It is given the symbol $J_\mathrm{1s2s}$ (the subscripts indicating the orbitals that the electrons are in). Because this is an expectation value for a potential energy arising from repulsion, it is necessarily a positive quantity.

The second and third terms are actually also the same, again because the variables are just dummy variables and you can interchange the labels. This integral, however, does not have a physical interpretation; it is a purely quantum mechanical effect. It is called the exchange integral, and is given the symbol $K_\mathrm{1s2s}$. It's a bit harder to show, but the exchange integral is also a positive quantity (there is a question on it here). However, its magnitude is not as large as $J_\mathrm{1s2s}$.

Putting it all together (leaving out subscripts for clarity) we have

$$\left\langle\Psi_\mathrm{triplet}\middle|\frac{1}{r_{12}}\middle|\Psi_\mathrm{triplet}\right\rangle = \frac{1}{2}(J-K-K+J) = J-K$$

and if you were to perform the same calculation on the singlet spin wavefunction (not very difficult, just replace the minus signs with a plus!) you would get 

$$\left\langle\Psi_\mathrm{singlet}\middle|\frac{1}{r_{12}}\middle|\Psi_\mathrm{singlet}\right\rangle = \frac{1}{2}(J+K+K+J) = J+K$$

Since $J &gt; K &gt; 0$, the expectation value of the energy for the singlet state is higher.

But I have more than two electrons. What now?

Well, you just have to do a ton more ugly maths and perhaps a few hundred more integrals. But the idea is the same. For every pair of electrons with parallel spins, there is a corresponding exchange integral that reduces the expectation value of the energy.

Therefore, in a $\mathrm{3d^5}$ manganese atom (we ignore the closed shells) there are

$$\frac{5!}{2!(5-2)!} = 10$$

different pairs of d-d electrons with parallel spins. Therefore if you evaluate the expectation value of the energy, there will be a $-10K_\mathrm{3d3d}$ term, which is stabilising. This is precisely the exchange energy effect that Geoff describes in his answer.