Why is laser light described by a coherent state? This is a follow-up to this recent answer by Wouter to this related question from 2015, and a comment by Emilio Pisanty underneath.
I have read the papers by Mølmer, Bartlett et al., Wiseman, and parts of van Enk and Fuchs referenced therein. I have a clear idea of the justification (even derivation) of the mixed state
$$\frac1{2π} \int |α e^{iφ}\rangle \langle α e^{iφ}| \mathrm{d}φ = e^{-μ} \sum_n \frac{μ^n}{n!} |n\rangle \langle n|.$$
But I want to agree with Bartlett et al.'s "factist" and see the emergence of the pure state
$$|\sqrtμ\rangle = e^{-μ/2} \sum_n \sqrt{\frac{μ^n}{n!}} |n\rangle.$$
Here's what Wiseman has to say about it:

Thus, the phase reference laser beam in the teleportation experiment is in a coherent state.
It is not a coherent state of unknown phase; it is the coherent state $|μ\rangle$, with zero phase (relative to itself).

and

In particular, if the laser itself is the clock then by definition it is coherent with respect to itself so it should be described by a pure state $|μ\rangle$ of zero phase.

I don't understand these claims. What does the author mean by "then by definition it is coherent with respect to itself"? Surely within coherence theory (the $g^{(n)}$ functions) an ideal laser is coherent to all degrees but I am quite sure that is a property of the Poissonian mixture as well, which would not prioritize the pure state in any way.
Furthermore, the arguments of virtually everyone are refuted in the Bartlett et al. paper, including Wiseman's quote above and van Eck and Fuchs' reference to the de Finetti theorem, but I can not infer the particular form of the pure state from any part of the talk about the reference frames, only that the two lead to the same measurable results. (I am very sure I am just lacking some context there.) I don't need to show that the coherent state works, or that it is not necessary for explaining the experiments. I want to derive it.
My question is simple: how, of all possible states which share the same photon statistics and coherence degrees, do I extract exactly the formula for the coherent state?
The discussion has grown just far too long and deep to follow in the original literature, so I'm asking for a quick summary from someone who's already familiar with it.
 A:  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.

A: After hours and hours of pondering this problem, here's an idea what taking the laser state "with respect to itself" might have meant.
Suppose that rather than trying to derive, one wants to measure what the state of the field is. One has the usual options of photon-number-resolving detection for $⟨n|ρ|n⟩$, homodyne detection for $⟨x_φ|ρ|x_φ⟩$, or eight-port homodyne detection for $⟨α|ρ|α⟩$. We already know the result of the first one, but the other two require a reference wave to compare with. We don't have it, so the best we can do is to use the very same laser instead, basically to measure it using itself.
Note that this is inherently circular logic. Homodyne detection and friends are known to work if a strong coherent state is in the local oscillator input. We substitute it with a laser in order to see... that it is in a strong coherent state! I guess that's what is meant by the "coherence with respect to itself", and ultimately by the "convenient fiction" of Mølmer: in fact we don't prove that the state is a coherent state, but using it as a reference we establish a basis in which it happens to be.
Nevertheless, carrying out this idea this is what we get: Consider a balanced homodyne detector where both the signal and the local oscillator are created by splitting the same input beam which is in the state
$$ρ = e^{-μ} \sum_{n=0}^∞ \frac{μ^n}{n!} |n⟩⟨n|.$$
This basically creates a Mach–Zehnder interferometer, the only difference being that a difference in photocounts is measured rather than the counts in one output alone. Elementary algebra shows that a Fock state $|n⟩$ produces difference photon number statistics
$$P(n-2k|n) = \binom{n}{k} \sin^{2k} φ \cos^{2(n-k)} φ$$
for the homodyne relative phase $φ$, and thus
$$\begin{aligned}P(m) &= \sum_{n=m}^∞ P(m|n) P(n) = \sum_{k=0}^∞ P(m|m+2k) P(m+2k) \\ &= \sum_{k=0}^∞ e^{-μ} \frac{μ^{m+2k}}{(m+2k)!} \binom{m+2k}{k} \sin^{2k} φ \cos^{2(m+k)} φ\end{aligned}$$
if the state is the mixture $ρ$.
Now, assuming that our homodyne works in this regime as intended and provides access to the distribution of the quadrature $x_φ$, it can be shown with little math and some approximations for large $μ$ that for $x_φ := m/\sqrt{2μ}$ the distribution becomes similar to a continuous normal distribution with mean $\sqrt{μ/2} \cos φ$ and standard deviation $1/\sqrt2$. (I may have some 2's the wrong way, but who cares.)
This is the result which would be measured for the pure state $|μ⟩$, and that $\tilde ρ = |μ⟩⟨μ|$ is the unique state consistent with this measurement statistics can be argued using the inverse Radon transform.

In the light of the above references, although I admit that I did not re-read them, perhaps the idea above (of ascribing to the field a state obtained by an unjustified measurement method) is not as perverted as it may feel.
Note that while measurements in the Fock basis are unambiguous, the relative phases of the basis vectors are our free choice. So the only question would be whether the “actual” state is a superposition of all of them, rather than a mixture – the phases could then be coerced to any pattern, including all zeroes, solely by specification of the basis.
It could be said that all we do in the above measurement scheme is fix these phases by measuring relative to the laser – the only, although significant, controversial point is calling it homodyne detection, and using its results in subsequent calculations as such. These oscillations may, and will, fluctuate in time in result of finite coherence time of the source, making the basis time-dependent, or rather, stochastic.
One could take a step outside, sample the basis vectors once and then just evolve them by the free Hamiltonian. This would create a different basis, in which the $P$ distribution would diffuse in time by Fokker-Planck to ultimately reach the stationary form of the well-known annulus of radius $\sqrtμ$ around zero, which describes the mixture $ρ$.
Indeed, I have become pretty sure that this is exactly what the most recent of the listed papers meant by the internal and external reference frames.
