Well, theory is of course consistent and you get an eigenstate with eigenvalue 0. However, as this is a zero ket, it will have exactly 0 weight in any experiment and whatever follow-up operation (such as photon creation) you perform, all of this processes will have 0 weight.
I just want to point out that photon number measurements usually do not implement the photon number operator. Instead one rather applies an operator that does something like (\hat{c}^\dagger a)^n. You destroy all of the n photons you have and instead create different excitations - usually photoelectrons or something like that, so strictly speaking the creation operator \hat{c}^\dagger will create a photoelectron in a slightly different mode for every photon. This is obviously not an implementation of the photon number operator as the photons are gone afterwards.
If one wants to implement the photon number operator, this is nontrivial and one has to do it in the manner as outlined in the manuscript above. However, from the experimental point of view, it is hard to implement the annihilation operator. You want to annihilate exactly one photon, not 2 photons or all of them, but exactly one. The usual way - as in the manuscript I linked - to implement the operator is a beam splitter, usually a glass plate, that splits off only a tiny fraction of the beam followed by a single-photon sensitive detector that verifies the absorption process. If the photon number in the original field is low enough, the probability that two or more photons were reflected by the beam splitter simultaneously may be made negligible.
Assuming for simplicity that the state could be described by a diagonal density matrix before, where the probability to have n photons was given by p_n, the probabilities to find n' photons in the photon-subtracted light field will be given by p_{n'}=p_n|\langle n'|\hat{a}|n\rangle|^2. This obviously only yields contributions for n'=n-1 and you get p_{n'=n-1}=p_n n. So, the weight of finding n photons in the photon-subtracted state is given by the probability that you had n+1 photons in the initial state times the number of photons. This makes sense as the operator you applied is the single photon annihilation operator and when placing a single beam splitter even if you had the same initial probability to find e.g. 1 or 5 photons inside your beam, it will happen 5 times as often that a single photon gets reflected from the n=5 state as compared to the n=1 state. Thus, the detector will also click 5 times as often. This rescaling of probabilities is at the core of many non-intuitive effects in quantum optics, such as doubling of the photon number inside a thermal beam if you subtract a photon. However, it also means that the vacuum contribution of the initial state will not contribute anything here. The probability that one photon gets reflected from the beam splitter resulting in a detector click is of course exactly 0.
This was the implementation of the annihilation operator. Implementing the creation operator now would be "simple". One would just have to add a single photon source that adds a photon to be the beam as soon as the detector clicked. In the manuscript I linked before this was done in a stochastic manner using heralding an an entangled photon pair. When this is done, you of course get the correct photon numbers out of it. However, in terms of finding eigenstates in a measurement and therefore in terms of state preparation, the behavior is somewhat different. If I consider a standard experiment such as a spin measurement or a position measurement on some mixed state, I expect the different eigenvalues to show up with a relative probability given by their weights. This is not as easy for the photon number operator. Here, for example, the probability to find the system in the vacuum state after applying the photon number operator is a constant 0. This is of course fully consistent, but has important consequences for experiments and e.g. preparation of states.
