This is a very easy question: I'm in need of some intuition on the fact that, e.g. thermal sources, produce bunched photons. It is very easy to "undertand", without any quantum mechanics, why single emitters produce antibunching. Is there any such explanation for bunching?
2 Answers
Ilya's answer describes "photon bunching" or more precisely intensity fluctuations in a laser or a medium with amplified spontaneous emission. In a real thermal source the situation is completely different. In contrast to the described stimulated emission processes the emission from two single photon emitters in a thermal light source is uncorrelated both in time and in phase. This yielding bunched photons seems counterintuitive if one imagines photons as bullets flying towards your detector on a classical trajectory. But if you consider the phase of the photons and the superposition principle in electrodynamics, then you see that it's exactly the fact that the phases of individually emitted photons are uncorrelated which leads to thermal photon statistics and hence, bunching.
Take the thermal light statistics expressed in terms of the photon creation and annihilation operators $\hat{a}^{\dagger}$, $\hat{a}$: $$ \left\langle \left( \hat{a}^{\dagger} \right)^l \left( \hat{a} \right)^k \right\rangle_{\text{therm}} = \delta_{kl} k! \bar{n}^k $$ $\bar{n}$ is the mean number of photons per mode. Now if the photon creation and annihilation operators $\hat{a}^{\dagger}$, $\hat{a}$ are expressed as superposition of the creation and annihilation operators of $N$ photons from individual atoms $\hat{a}_i^{\dagger}$, $\hat{a}_i$ with uncorrelated phases $\phi_i$ $$ \hat{a}^{\dagger} = \frac{1}{\sqrt{N}} \sum_{i=1}^{N} \hat{a}_i^{\dagger} e^{-i \phi_i} \\ \hat{a} = \frac{1}{\sqrt{N}} \sum_{i=1}^{N} \hat{a}_i e^{i \phi_i} $$ you get $$ \left\langle \left( \hat{a}^{\dagger} \right)^l \left( \hat{a} \right)^k \right\rangle = \frac{1}{\sqrt{N^{k+l}}} \sum_{i_1, ..., i_l = 1}^{N} \sum_{j_1, ..., j_k = 1}^{N} \left\langle \hat{a}_{i_1}^{\dagger} e^{-i \phi_{i_1}} \cdot ... \cdot \hat{a}_{i_l}^{\dagger} e^{-i \phi_{i_l}} \cdot \hat{a}_{j_k} e^{i \phi_{j_k}} \cdot ... \cdot \hat{a}_{j_1} e^{i \phi_{j_1}} \right\rangle . $$ Calculating the expectation values by averaging over many realizations of phases $\phi_i$ many summands go to $0$ because $\left\langle e^{i \phi} \right\rangle_{\phi} = 0$. The only summands that survive the averaging are those, where all phase terms cancel out, i. e. if for every $e^{\phi_i}$ there is a conjugate phase $e^{-\phi_i}$. This can obviously only be the case for $l = k$. This is where the Kronecker Delta $\delta_{kl}$ comes from. As the phases of photons from individual atoms are uncorrelated the conjugate phases can only originate from the creation $\hat{a}_i^{\dagger} e^{-i \phi_i}$ and annihilation term $\hat{a}_i e^{i \phi_i}$ corresponding to the same atom.
As all nonzero expectation values in the sum contribute equally (namely with $\left\langle \hat{a}_i^{\dagger} \hat{a}_i \right\rangle^k = \bar{n}^k$) we just need to find the number of nonzero expectation values. This is given by the number of possibilities to map the $i_a$s onto the $j_b$s ($a$, $b$ $\in \left\{ 1, ..., k \right\}$) bijectively. Multiple occurances of the same index are not possible, because each single photon emitter can only contribute one photon in each realization.
With $\frac{N!}{(N-k)!} k!$ being the number of nonzero summands, we arrive at $$ \left\langle \left( \hat{a}^{\dagger} \right)^k \left( \hat{a} \right)^k \right\rangle = \frac{1}{\sqrt{N^{2k}}} \frac{N!}{(N-k)!} k! \left\langle \hat{a}_i^{\dagger} \hat{a}_i \right\rangle^k = \frac{1}{N^k} \frac{N!}{(N-k)!} k! \bar{n}^k . $$ Combining this result with the previously obtained Kronecker Delta and taking the limit of extremely many single photon emitter $N$ it yields $$ \lim_{N \to \infty} \delta_{kl} \frac{1}{N^k} \frac{N!}{(N-k)!} k! \bar{n}^k = \delta_{kl} k! \bar{n}^k , $$ which is exaclty the thermal light statistics.
As you can see the thermal light statistics (including bunching in all orders $k$) arise from the assumption of uncorrelatedly emitting single photon sources. As an intuitive picture of these intensity fluctuations you can imagine that in some realizations (or points in spacetime) the photons' phases interfere contructively so that you get high intensity and in other realizations they interfere destructively and you get less intensity.
This kind of chaotic interference can also be observed in macroscopic systems. Consider a laser beam impinging on a phase mask (e. g. piece of ground glass): Then at each position in the laser spot the light accumulates a different phase $\phi_i$. After the mask the parts of the beam interfere and form a speckle pattern. This speckle pattern corresponds to one realization of multiple sources emitting with individual phases $\phi_i$. If you change the phase mask, for example by moving it such that the laser hits a different spot on the ground glass you get a different realization of phases $\phi_i$ and hence a different speckle pattern. Imagine moving the ground glass continuously – this would result in a continuously changing speckle pattern. A point-like detector will then sometimes be sitting on a brighter spot and at other times in a dimmer region.
In real thermal light this is somewhat similar. The difference is that the interference pattern changes much faster, so that your eye cannot resolve the speckles anymore. And also the average number of photons per mode $\bar{n}$ ($\approx$ number of photons in each speckle spot) is usually lower than 1, so that even a high-speed camera would not see a speckle pattern but instead single photons arriving at random positions and random times. Only if you then calculate correlation functions of the photon arrival spacetime coordinates you will find that the probability of detecting a second photon on the same spot as a directly preceeding photon is increased by a factor of $2$ compared to the probability of detecting a photon on any other spot on average.
Here is probably the simplest argument that I have heard of. Detection of a photon from a thermal source gives a rise to a probability to detect several more in a short interval of time due to stimulated emission. Assume that you have some atoms in a medium that emits light, and they are in an excited state. If you know that one atom emitted a photon, this will increase the probability to emit a photon for its neighbors. So you will have a bunch of photons leaving this body and many atoms in a ground state, that need some time to "recharge". While in a laser you create some artificial conditions when there are always "enough" of excited atoms, moreover there is a feedback: field goes stronger -> lower level of population inversion -> field weakens and vice versa. Therefore there is no bunching/antibunching for laser.