Why does classical light always result in super-Poissonian statistics? It is a well-known result that classical light (which I take here to mean mixtures of coherent states) cannot produce sub-Poissonian photon-counting statistics, with a single beam of coherent light corresponding to a Poissonian photon-counting statistics (as discussed for example here), and other kinds of non-quantum light corresponding to super-Poissonian statistics.
However, I have never seen this fact proven formally. Usually, texts show how some common kinds of classical light, such as thermal light, result in super-Poissonian statistics, and how quantum states can produce sub-Poissonian ones, but they do not tackle the general case.
More specifically, consider a state which is a mixture of coherent states. This corresponds to a photon counting probability $P(n)$ of the form
$$P(n)=\sum_\lambda p_\lambda P_\lambda(n),$$
with $\sum_\lambda p_\lambda =1$, and $P_\lambda(n)$ being the Poisson distribution with expected value $\lambda$:
$$P_\lambda(n)\equiv e^{-\lambda}\frac{\lambda^n}{n!}.$$
A super-Poissonian distribution is characterised by the property that the variance is greater than the expected value, that is, $\sigma^2\ge\mu$. More precisely, in the considered case this means
$$\sum_n(n-\mu)^2P(n)\ge \mu,\quad \mu\equiv\sum_n nP(n).$$
Can this property be shown in full generality, without making reference to specific types of light?
 A: Let us compute the first moments of $P(n)$:
$$\mu\equiv\sum_n nP(n)=\sum_n n\sum_\lambda p_\lambda P_\lambda(n)=\sum_\lambda p_\lambda \lambda,$$
where I used the property of the Poisson distribution $\sum_n n P_\lambda(n)=\lambda$.
Similarly, we have
$$\sum_n n^2 P(n)=\sum_\lambda p_\lambda \lambda(\lambda+1),$$
where I used $\sum_n n^2 P_\lambda(n)=\lambda(\lambda+1)$.
The variance $\sigma^2$ of the distribution thus reads
$$\sigma^2\equiv\sum_n (n-\mu)^2 P(n)=\sum_\lambda p_\lambda \lambda(\lambda+1)-\mu^2,$$
and finally the difference between variance and expected value, $\sigma^2-\mu$, is
$$\sigma^2-\mu=\sum_\lambda p_\lambda\lambda(\lambda+1)-\mu(\mu+1).\tag1$$
Defining $f(\lambda)\equiv\lambda(\lambda+1)$, (1) can be written as
$$\sigma^2-\mu=\sum_\lambda p_\lambda f(\lambda)-f\Big(\underbrace{\sum_\lambda p_\lambda \lambda}_{\mu}\Big).$$
The conclusion $\sigma^2-\mu\ge0$ now follows from $f$ being convex, together with Jensen's inequality.
This proves that an arbitrary mixture (convex combination) of Poissonians gives a super-Poissonian distribution satisfying $\sigma^2\ge\mu$.
