A true but misleading identity
Consider a single mode, and let $|n\rangle$ be the state with $n$ photons in that mode. A coherent state has the form
$$
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\da}{a^\dagger}
|re^{i\phi}\ra
=
\sum_{n\geq 0} c_n (re^{i\phi})^n|n\ra
\tag{1}
$$
where the coefficients $c_n$ depend only on $n$, not on the real numbers $r$ or $\phi$. Use
$$
\frac{1}{2\pi}\int_0^{2\pi} e^{i(n-m)\phi}= \delta_{n,m}
$$
to get the identity
$$
\frac{1}{2\pi}
\int_0^{2\pi}d\phi\ |re^{i\phi}\ra\,\la r e^{i\phi}|
=
\sum_{n\geq 0} |c_n|^2 r^{2n}|n\ra\,\la n|.
\tag{2}
$$
This identity shows that the density matrix representing a coherent state with known amplitude $r$ but unknown phase is the same as the density matrix representing an unknown number eigenstate with a particular probability distribution for the number $n$.
That's misleading, though, because the density matrix for the light by itself does not tell the whole story. It does not convey everything we know about how the light was produced. That additional information can have observable consequences.
Analogy
To see why (2) does not tell the whole story, consider a simpler analogy.
Suppose we measure the spin of an electron along a particular axis. The possible outcomes are $|u\ra$ and $|d\ra$, for "up" and "down". But suppose that we got distracted and didn't notice which outcome occurred. We know that the measurement occurred, and we know what the axis was, so we know what the two possible outcomes are, but we don't know which of those two outcomes was actually obtained.
What density matrix should we use to represent the electron's spin-state if we want to make predictions about subsequent spin measurements? The answer to that question is
$$
\rho\propto |u\ra\,\la u|+|d\ra\,\la d|,
\tag{3}
$$
which is proportional to the identity matrix. That's the same density matrix that we should use no matter what the first measurement's axis was, a situation that is analogous to the identity (2). This is as it should be, because if we don't know the outcome of that first measurement, then we can't predict anything about the outcome of the next measurement of the electron's spin by itself.
However, the density matrix for the electron's spin does not tell the whole story. We have more information, because we know what measurement was done, even though we don't know the outcome. To account for this, we can use a density matrix for the combined system that includes the electron's spin and the device that measured it, like this:
$$
\rho'\propto \sum_M
\big(|u, M_u\ra + |d,M_d\ra\big)
\big(\la u, M_u| + \la d,M_d|\big)
\tag{4}
$$
where $M_{u/d}$ represents the state of the measuring device in the wake of the outcomes $u/d$ and the sum over $M$ is a sum over unknown microscopic details of the measuring device prior to the measurement. Calling it a measurement implies $\la M_u|M_d\ra=0$, so taking a partial trace over $M$ gives (3). But unlike the electron-only density matrix (3), the density matrix (4) for the combined system does depend on the measurement axis.
That additional information doesn't help us predict the outcome of a subsequent measurement of the electron's spin by itself, but it might help us predict other things. As an example, suppose that the process of measuring the electron's spin produces a photon whose polarization is $V$ or $H$ according to whether the spin-measurement outcome was $u$ or $d$. The reduced density matrix (3) doesn't know anything about this correlation, but the more complete density matrix (4) does, and that correlation obviously has observable consequences.
Laser light
We could derive something similar to (4) for the light from a laser, where now $M_\circ$ is something about the device that produced that light, but that would require keeping track of microscopic details of the device. We can do something easier instead. We can consider the density matrix of the light by itself, but instead of considering only the final state, we can consider how the state evolves in time under the influence of the device that is producing the light. Then we can ask which initially-pure density matrices remain most nearly pure for the longest time.
To see why that works, consider the spin-measurement analogy again. The initially pure density matrices $|u\ra\,\la u|$ and $|d\ra\,\la d|$ remain pure when the spin is measured along the $u/d$ axis, but any other initially pure density matrix becomes mixed as soon as the spin is measured along the $u/d$ axis. In this way, the information about the measurement-axis is revealed by the time-dependence of the spin-only density matrix.
Similarly, information about the device that produced the laser light is implicit in the time-dependence of the light-only density matrix: the states that would appear in the combined-system density matrix (analogous to (4)) should correspond to the pure states that remain most nearly pure for the longest time. The details are worked out in
- Gea-Banacloche (1998) "Emergence of classical radiation fields through decoherence in the Scully-Lamb laser model," Foundations of Physics 28: 531-548.
The conclusion is previewed on page 533:
Near threshold, the state [that remains most nearly pure for the longest time] is very nearly a coherent state, whereas high above threshold it is found to be slightly squeezed in intensity, i.e., with slightly sub-Poissonian photon statistics.
...and restated later on page 546:
near threshold the longest-lived state of the field in a laser is very nearly a coherent state, and well above threshold a sub-Poissonian (i.e., nonclassical) state — although still "quasiclassical" in the sense of having a fairly well-defined phase and amplitude.
Page 546 also acknowledges a limitation of the analysis:
The calculations [above] have concerned themselves with the initial decay rate of various possible field states, without actually attempting to track how that decay rate may, in fact, change with time. In their related article on decoherence in a linear harmonic oscillator [ref], Zurek and co-workers find at least one instance in which squeezed states show a smaller initial decay rate than coherent states, yet after a finite time elapses the coherent states are found to retain more of their initial purity than the squeezed states. The characteristic time scale for this crossover to take place is the oscillator period. ... [The] main difference... is not so much between coherent and/or slightly squeezed states, but between quasiclassical states, ...and the very nonclassical energy eigenstates (number states) which decay at the much faster rate...
I won't repeat the calculation, but here are some highlights. In the Scully-Lamb model, the evolution of the light-only density matrix is given by
\begin{align}
\frac{d\rho}{dt}
=&\ -\frac{A}{2}(a\da\rho+\rho a\da - 2\da\rho a)
\\
&\ +\frac{B}{8}(a\da a\da\rho + 6a\da\rho a\da
+\rho a\da a\da - 4\da\rho a\da a- 4\da a\da\rho a)
\tag{5}
\\
&\ -\frac{C}{2}(\da a\rho + \rho\da a-2 a\rho\da)
\end{align}
where $a,\da$ are the photon annihilation/creation operators for a single mode, the term with coefficient $A$ accounts for stimulated and spontaneous emission, $B$ is a saturation coefficient, and $C$ is the cavity decay rate. Near threshold, $A\sim C$. The paper calculats $\la n|\rho(t)|n\ra$ for a number eigenstate and finds that, near threshold, the "lifetime" of such a state is $\sim 1/(nC)$. The paper also calculates $\la\alpha|\rho(t)|\alpha\ra$ for a coherent state with the same amplitude (same mean number of photons) and finds that the lifetime is $\sim 1/C$, a factor of $n$ longer than the lifetime of the number eigenstate. Since the mean value of $n$ is very large for a bright laser, this is a very significant difference in lifetimes.
This result can be understood intuitively. The process of producing laser light involves emission and absorption. Emission and absorption of individual photons is very disruptive to a number-eigenstate, but it hardly has any effect on a coherent state. More precisely: adding or removing one photon to/from a number eigenstate makes the new state orthogonal to the original, but adding or removing one photon to/from a large-amplitude coherent state hardly changes the state at all: the new state is approximately proportional to the original. (If a photon is removed, the new state is exactly proportional to the original. This is the defining property of coherent states.)
This leads us to expect that something similar to coherent states should be the longest-lived states, as confirmed by the calculation in the paper. Thus something similar to coherent states correspond to the states $u/d$ in the spin-measurement analogy (4).
A simplified version
To supplement the intuition given above, consider this simplified version of (5), in which only the cavity-decay terms are retained:
$$
\frac{d\rho}{dt}
= -\frac{C}{2}(\da a\rho + \rho\da a-2 a\rho\da).
\tag{6}
$$
Now the calculations are easy enough to do in your head. For a number eigenstate $\rho=|n\ra\,\la n|$ such that $\da a|n\ra=n|n\ra$, the decay rate is
$$
\frac{\la n|\dot\rho|n\ra}{\la n|\rho|n\ra} = -nC,
\tag{7}
$$
and for a coherent state $\rho=|\alpha\ra\,\la\alpha|$ such that $a|\alpha\ra=\alpha|\alpha\ra$, the decay rate is
$$
\frac{\la \alpha|\dot\rho|\alpha\ra}{\la \alpha|\rho|\alpha\ra} = 0.
\tag{8}
$$
In this oversimplified model, number eigenstates lose their purity at a rate proportional to $n$, but a coherent state tends to be more robust, at least initially.
I'll finish with one more excerpt from the paper:
it is, I think, a remarkable result that the very particle-like processes of emission and absorption of individual photons end up selecting for the kinds of states with the most wave-like properties — the coherent, or near-coherent, states.