The only correct way known is QED, and there is no photon wave function in the usual sense, because a photon has no position operator to begin with. Also photons are easily destroyed and created in interactions. So photon number is not conserved. Another point to consider is gauge invariance.
The most pragamatic approach is to start with classical Maxwell electrodynamics, introducing the potentials and fix the gauge completely by, e.g., choosing the Coulomb gauge.
To define "photons" you consider the free em-field case first and canonically quantize it. The Coulomb gauge in this case leads to the radiation gauge, where the classical (!) potentials are subject to the constraints
$$\vec{\nabla} \cdot \vec{A}=0, \quad \Phi \equiv A^0=0.$$
Now you quantize the fields canonically with the qualification that the canonical equal-time commutation relations must be made consistent with the transversality condition.
Finally this leads to the mode decomposition of the quantized field ##\vec{A}## with only two physical field-degrees of freedom. The most convenient choice of a single-photon basis is the momentum-helicity eigenbasis with the corresponding annihilation and creation operators ##\hat{a}(\vec{k},\lambda)## and ##\hat{a}^{\dagger}(\vec{k},\lambda)## with ##\vec{k} \in \mathbb{R}^3## and ##\lambda \in \{-1,1\}##. These operators follow the usual commutation rules (for bosons, which is necessary to have a local (microcausal) QFT with the Hamiltonian bounded from below),
$$[\hat{a}(\vec{k},\lambda),\hat{a}^{\dagger}(\vec{k}',\lambda')]=(2 \pi)^3 \delta^{(3)}(\vec{k}-\vec{k}') \delta_{\lambda \lambda'}.$$
A complete basis is the Fock basis, which is the (generalized) eigenbasis of the number operators ##\hat{N}(vec{k},\lambda)=\hat{a}^{\dagger}(\vec{k},\lambda) \hat{a}(\vec{k},\lambda)##.
To calculate a photon-number distribution you need to specify the state. For thermal photons, described by the canonical ensemble, the state is given by
$$\hat{\rho}=\frac{1}{Z} \exp \left (-\frac{\hat{H}}{T} \right), \quad Z=\mathrm{Tr} \exp \left (-\frac{\hat{H}}{T} \right).$$
The Hamiltonian is given by
$$\hat{H}=\sum_{\lambda \in \{-1,1\}} \int_{\mathbb{R}^3} \mathrm{d}^3 k |\vec{k}| \hat{N}(\vec{k},\lambda).$$
The photon-number distribution then turns out to be the Planck distribution (as it should be), i.e.,
$$f=\frac{2}{\exp \left (\frac{|\vec{k}|}{T} \right)-1}.$$
The ##2## comes from the two helicities.
Another example is a coherent state, which you get when assuming the source being (approximately) described by a classical charge-current distribution ##(\rho,\vec{j})## (the socalled "hemiclassical approximation", to be distinguished from the "semiclassical approximations, where the charges are quantized but the em. field kept classical). That's also how you describe quantum-field theoretically the emission of an antenna (e.g., the usual Hertzian dipole approximation treated in the standard classical-electrodynamics lecture). Considering harmonic time dependence you get a Poisson distribution for the photon number.
Last but not least there are also Fock states of definite photon numbers, which you can prepare by simply using an excited atom and wait till it deexcites and emits precisely one photon in a dipole transition or, more efficiently, by parametric downconversion using a laser (emitting never photon Fock states but coherent states) and a BBO crystal to produce entangled two-photon Fock states.
