In quantum mechanics, a coherent state of a quantum harmonic oscillator (QHO) is an eigenstate of the lowering operator. Expanding in the number basis, we find that the number of photons in a coherent state follows a Poisson distribution.
Is there a simple and intuitive reason why this fact holds?
$\newcommand{\ket}[1]{\lvert #1 \rangle}\newcommand{\Ket}[1]{\left| #1 \right>} % $Because you can use beamsplitters to split a coherent state into a tensorial product of many independent low photon number coherent states.
If you send $\ket{\alpha}$ on a beamsplitter of transmission coefficient $t$ and reflection coefficient $r$ (with $|r|^2+|t|^2=1$), you obtain the product of two independent coherent states $\ket{t\alpha}\otimes\ket{r\alpha}$. This property characterizes coherent states, since any other input state leads to entanglement in the output of the beamsplitter.
Since the output state is a product state, the statistics of any measurement done at one output is independent from the ones of a measurement performed at the other output. Furthermore, since the beamsplitter is a passive component, the total number of photons of the input state $\ket{\alpha}$ is the sum of number of photons at the outputs.
Now, you can also add beamsplitters at the outputs, and construct a tree of beam splitters, with $N\gg|\alpha|^2$ balanced outputs, transforming the input coherent state $\ket{α}$ into the product of $N$ coherent states $\Ket{\tfrac{\alpha}{\sqrt{N}}}^{\otimes N}$. As before, the total number of photon is conserved, thus the statistics of the number of photons of $\ket α$ is the sum of the $N$ independent outputs, each having a small average photon number $\tfrac{|\alpha|^2}{N}$. When $N \to \infty$, the only distribution having this property is the Poisson distribution. QED.
Note that, in the reasoning above, the beamsplitters do not need to be actual object splitting beams. Anything that changes the basis of space-time modes does the job. In particular, let your coherent state be in the mode corresponding to a pulse of light. You can also “slice” the pulse into $N$ short time slices. This description is exactly equivalent to the beamsplitter above, and corresponds to the intuition formulated by @AccidentalFourierTransform and @ThomasS above about the independence of successive photon detection events.
In all the descriptions above, I have implicitly assumed that the other port of each beamsplitter is empty, that is receives the vacuum state $\ket0$. This crucial assumption is still present above when I “slice” the coherent state in many timeslices, the initial $N-1$ vacua being in spacetime modes, which are orthogonal to the original lightpulse.
While the accepted answer already nicely answers the question, I believe it can be nice to see more explicitly how exactly it is that we get the coefficients (and thus the Poissonian statistics) of a coherent state $\lvert\alpha\rangle$ from the sole requirement that, after a unitary evolution $\mathcal U$, the output state $\mathcal U\lvert\alpha\rangle$ is factorised over the different modes: $$\mathcal U\lvert\alpha\rangle=\bigotimes_k\lvert\psi_k\rangle.$$ Let us consider a generic initial single-mode state $\lvert\psi\rangle$, with coefficients $$\lvert\psi\rangle=\sum_{k=0}^\infty \frac{c_k}{\sqrt{k!}}\lvert k\rangle = \sum_{k=0}^\infty \frac{c_k}{k!}a^{\dagger k}\lvert 0\rangle.$$ Note that the choice of $c_k/\sqrt{k!}$ as coefficients (instead of a simpler $c_k$) is purely conventional (but it will turn out to ease the calculations afterwards). Let us consider a unitary evolution which splits the state into two different modes as follows $$a^\dagger\to r a_1^\dagger + t a_2^\dagger.$$ We could consider the more general case of a splitting into $N>2$ modes, but this will turn out to not be necessary.
Upon this unitary evolution, the state $\lvert\psi\rangle$ evolves into: $$\lvert\psi\rangle\to\sum_{k=0}^\infty \frac{c_k}{k!}(ra_1^\dagger+ta_2^\dagger)^k\lvert0\rangle = \sum_{k=0}^\infty c_k\sum_{l=0}^k\frac{(ra_1^\dagger)^l}{l!}\frac{(ta_2^\dagger)^{k-l}}{(k-l)!}\lvert0\rangle,\tag1$$ where we used the binomial formula to expand the $k$-th power of a sum of two terms. Let us now ask what are the terms in this sum that contain the $n$-th power of $a_1^\dagger$. The answer is readily seen to be $$\frac{(ra_1^\dagger)^n}{n!}\sum_{k=n}^\infty c_k\frac{(ta_2^\dagger)^{k-n}}{(k-n)!}=\frac{(ra_1^\dagger)^n}{n!}\sum_{m=0}^\infty c_{m+n}\frac{(ta_2^\dagger)^m}{m!}.$$ In other words, if we rearrange the terms of (1) to make the coefficients of the powers of $a_1^\dagger$ explicit, we can rewrite the final state as following (1) $$\lvert\psi\rangle\to\mathcal U\lvert\psi\rangle= \sum_{n=0}^\infty \frac{(ra_1^\dagger)^n}{n!} \sum_{m=0}^\infty c_{n+m}\frac{(ta_2^\dagger)^m}{m!}\lvert0\rangle.$$ Remarkably, this tells us that the output state $\mathcal U\lvert\psi\rangle$ is separable if and only if the coefficients satisfy $c_{n+m}=c_n c_m$, that is, if $c_n=\alpha^n$ for some $\alpha\in\mathbb C$. As a bonus, we see that when this is the case, the output state remains a product of coherent states over the different modes.
(1) More coincisely, the gist of the argument is the following equality $$\sum_{s=0}^\infty c_s(a+b)^s=\sum_{n=0}^\infty \frac{a^n}{n!}\sum_{m=0}^\infty c_{n+m}(n+m)!\frac{b^m}{m!}.$$
1. From properties of destruction operator
So, first you have to accept that $a|n\rangle = \sqrt{n}|n-1\rangle$. This is relatively easy to see because the matrix element for absorption of a photon by a two-level system (atom which goes from ground to excited state) is proportional to $\langle n-1|a|n\rangle$ and this must be proportional to the square root of the number of photons in the light mode because the probability for absorption must be proportional to the light intensity. So you need something like $a|n\rangle = \sqrt{n}|n-1\rangle$ (ignoring a possible phase factor).
Then, when you expand the coherent state in number states, $|\alpha\rangle = \sum_n c_n |n\rangle$ and put this into $a|n\rangle = \alpha|n\rangle$, you see that you need $c_n\sqrt{n}=\alpha c_{n-1}$. The result when lowering $c_n|n\rangle$ with $a$ must be the same as a multiplication of $|n-1\rangle$ with $\alpha$. As a consequence, $c_n =(\alpha/\sqrt{n})c_{n-1}$ and you are finished. Iterating this $n$ times yields $c_n = (\alpha^n/\sqrt{n!})c_0$. Normalization gives the value of $c_0$ and then you have $\langle n|\alpha\rangle=c_n$. Now you square the whole stuff and get the Poisson distribution.
So the point is that for large $n$, $\alpha/\sqrt{n}$ will always be smaller than 1. This is why the Poisson distribution decreases in this case. For small $n$, the opposite holds and the Poisson distribution increases.
2. Coherent state in phase space
There's an alternative picture. You know that a single-mode field is like a harmonic oscillator where the mode's quadrature operators play the role of position and momentum of the HO. Now, a coherent state is a wave packet that oscillates in the parabolic potential without changing its shape. There is no dispersion for this wave packet, it coheres (this is where the name coherent state comes from). The energy eigenstates of the HO (which correspond to the number states of the field mode) are static, they don't move. So, to construct a coherent state, you need to use a superposition of number states. And the weighting of number states in the superposition is the square of the probabilities of the Poisson distribution.
This is also not an intuitive physical explanation but it shines a little bit more light on the problem.
3. Coherent state and independent emission events
Another possiblity to get a physical understanding is the independence of the "emission" events. From this, the Poisson distribution is easily understood. What I don't see is the connection between the coherent state $|\alpha\rangle$ and the concept of statistically independent emissions. I think it's even counterintuitive. In the laser, the induced emission events (together with the resonator) create the coherent state. The statistically independent spontaneous emission events disturb the coherent state (phase fluctuations in the laser).
Who can help?
The Poisson distribution is derived statistically from an input of random occurances,
it expresses the probability of a given number of events occurring in a fixed interval of time and/or space if these events occur with a known average rate and independently of the time since the last event.1 The Poisson distribution can also be used for the number of events in other specified intervals such as distance, area or volume.
The Poisson distribution is an appropriate model if the following assumptions are true.
K is the number of times an event occurs in an interval and K can take values 0, 1, 2, …
Check for photons.
The occurrence of one event does not affect the probability that a second event will occur. That is, events occur independently.
Check, there is no photon-photon interaction, only superposition.
The rate at which events occur is constant. The rate cannot be higher in some intervals and lower in other intervals.
Check?
Two events cannot occur at exactly the same instant.
Check, it is the heisenberg uncertainty here.
The probability of an event in an interval is proportional to the length of the interval.
Check?
If these conditions are true, then K is a Poisson random variable, and the distribution of K is a Poisson distribution.
I have question marks where I do not know what a QHO is.
If it checks, then that is the reason Poisson is used. There are two examples for photons in the occurrences list.
Physically the energy spent to create photons is much smaller than the transferred energy in collisions of charged particles. In the zeroth-order approximation one can neglect the influence of the lost energy on the particle scattering and then the particle behavior becomes known at all times. Often it is referred as to a "classical current" $j(t)$. On the other hand, the radiation equation in this approximation is, roughly speaking, linear in the known current $j$: $\square A =j$ and so are the corresponding Fourier components. To be exact, one has to solve the QED equations in this approximation, see formula (24.27) in Akhiezer-Berestetski textbook:
In this approximation the field solution for each harmonics is an eigenstate of the lowering operator $a_{\omega}$. It simply means two things: 1) emitted photons do not hinder/favorise each other while emitting, and 2) there is always enough energy to create any number of photons. In other words, emission of one photon does not influence emission of another one at the same or any other time. Their emission is random. Now statistics of randomly emitted number of photons comes into play and you get the Poisson distribution.
As soon as you fix or limit from above the emitted energy, the Poisson distribution gets spoiled, especially for high frequencies and high photon numbers.
