Q:

How to justify "correct Boltzmann counting" (dividing by $N!$) without resorting to quantum mechanics?

Let a box of volume $V$ contain $N$ monatomic gas particles, with total energy $E$. By considering the "volume" of phase space consistent with our macro-observables ($E$, $V$, and $N$), the multiplicity of this macrostate is (from Ref.1)
$$
\Omega = \frac{V^N}{h^{3N}}\frac{2\pi^{3N/2}(2mU)^{(3N-1)/2}}{\Gamma(3N/2)} \,.
$$
However, this expression is "incomplete" as it leads to a non-extensive entropy and hence Gibbs paradox. The fix for this is dividing by $N!$ (or "correct Boltzmann counting"), and the justification is usually: because particles are fundamentally indistinguishable in quantum mechanics.
I don't like the quantum mechanical justification, and Ref.2 sums up my feelings

If this reasoning were correct, then the Gibbs paradox would imply
  that the world is quantum mechanical.  It is difficult to believe that
  a thought experiment from classical physics could produce such a
  profound insight.

Furthermore, I intuitively accept that the probability of a macrostate is proportional to the volume of phase space it covers. Dividing by $N!$ would "unfairly" reduce the probability of macrostates with larger number of particles. If $E$, $V$, and $N$ were all doubled, then the compliant region of phase space would increase drastically, more than just getting squared (i.e. doubling in entropy).
On the other hand, if entropy is not an extensive quantity, then neither is chemical potential $\mu$, and using the earlier "uncorrected" multiplicity, an infinite particle reservoir at constant temperature and pressure will have non-finite chemical potential, this can't be correct. (Related to: How do you experimentally measure the chemical potential of a gas?)
According to Ref.2,

systematic division by $N!$ for all particles of the same kind in
  mutually accessible states has no empirical consequences.

and it goes on to say that this is because dividing by $N!$ in closed systems do not affect the relative probabilities of macrostates (as $N$ is constant).
My understanding didn't really "click" after reading Ref.2, so I suppose I'm looking for a more detailed (or "intuitive") argument as to why correct Boltzmann counting has no empirical consequences.

Schroeder, D. (1999). An Introduction to Thermal Physics. p.71.
Versteegh, M., Dieks, D. (2011). "The Gibbs Paradox and the distinguishability of identical particles".

This question is similar to Why is the partition function divided by $(h^{3N} N!)$?, but I'm looking for a non quantum mechanical explanation specifically.

A:

After thinking about this for a bit, this is what I think could be an explanation.

Let
$$
\Omega(V, N) = \frac{V^N}{h^{3N}}\frac{2\pi^{3N/2}(3Nmk_BT)^{(3N-1)/2}}{\Gamma(3N/2)}\Delta p
$$
for some constant $T$ and where $\Delta p$ is some uncertainty in momentum. This is the multiplicity for a monatomic ideal gas at some constant temperature without correct Boltzmann counting.
I will first argue that the uncorrected multiplicity is actually valid, and that the Gibbs paradox can be averted.
For later convenience, we simplify the expression for $\Omega(V, N)$.
\begin{align*}
    \Omega(V, N)
        &= \frac{V^N}{h^{3N}}\frac{2\pi^{3N/2}(3Nmk_BT)^{(3N-1)/2}}{\Gamma(3N/2)}\Delta p \\
        &= \frac{(3N/2)^{(3N-1)/2}}{\Gamma(3N/2)}\left[\frac{(2\pi mk_B T)^{3/2}}{h^3}V\right]^{N}\sqrt{\frac{2}{mk_BT}}\Delta p
\end{align*}
Let
\begin{align*}
    f(N) &= \frac{(3N/2)^{(3N-1)/2}}{\Gamma(3N/2)} \\
    \lambda &= \frac{h}{\sqrt{2\pi mk_B T}} \\
    \delta &= \sqrt{\frac{2}{mk_BT}}\Delta p
\end{align*}
then
\begin{align*}
    \Omega(V, N)
        &= f(N)\left(\frac{V}{\lambda^3}\right)^{N}\delta
\end{align*}
As per the usual setup for Gibbs paradox, consider the system below

                

where the left and right boxes are in thermal equilibrium but separated by a divider, and their volumes and number of particles are denoted by $V_1, N_1$, and $V_2, N_2$ respectively. We also constrain $N_1/V_1=N_2/V_2$, so both boxes have equal pressures.
It is then said that the total multiplicity of the system is
$$
\Omega^{(?)}_\text{sys} = \Omega(V_1, N_1)\Omega(V_2, N_2) \,,
$$
but I would argue that this is incorrect.
Suppose there is only one particle in each box, the system would then be characterised by the phase space coordinates $(\vec{q}_1,\vec{p}_1,\vec{q}_2,\vec{p}_2)$ where $\vec{q}_i$ is the position of particle $i$ and $\vec{p}_i$ is the momentum of particle $i$. The quantity $\Omega^{(?)}_\text{sys}$ would actually correspond to the multiplicity of a phase space distribution where the positions are distributed like so

                             

However, this represents only a subset of the states that are consistent with our macrovariables (which was, one particle in each box). The position distribution should actually be

                             

Because the particles are macroscopically identical (our macrovariables do not distinguish between them), there is no reason to believe that our system exists in only a subset of the available microstates as suggested by fig. 2. Of course, the system is definitely in one of the square regions in the position distribution, it cannot evolve between the two. If we wanted to insist that the system is in a specific square, we must have a way to prescribe this configuration using our macrovariables, but this is equivalent to requiring the particles be macroscopically distinguishable.
Another reason why $\Omega^{(?)}_\text{sys}$ is not the correct multiplicity is because it was calculated by multiplying the multiplicity of the left box subsystem, with that of the right box subsystem. This only works if the microstates of the total system can in some sense be factored into two uncorrelated "sub-microstates". If we look at fig. 3, the microstates cannot be factored into independent subparts. If $\vec{q}_1$ is a position in the left box, then $\vec{q}_2$ must be a position in the right box.
The system multiplicity should really be (in the one particle per box case)
$$
\Omega_\text{sys} = \Omega(V_1, 1)\Omega(V_2, 1) \times 2 \,,
$$
where we multiply by two to account for the correct volume of phase space. For a general number of particles per box, it is
$$
\Omega_\text{sys} = \Omega(V_1, N_1)\Omega(V_2, N_2) \times  {N_1+N_2 \choose N_1,N_2}\,,
$$
where we have instead multiplied by a multinomial coefficient. Intuitively, we are saying that if $(\vec{q}_i,\vec{p}_i)$ is in a region of phase space compliant with our macroscopic prescription, then so will any permutation of the position-momentum pairs.
We now remove the divider between the boxes and let the particles mix. The system multiplicity is now
$$
\Omega_\text{sys}' = \Omega(V_1+V_2, N_1 + N_2) \,.
$$
This turns out to be equal (kind of) to the multiplicity before the divider was removed $\Omega_\text{sys}$, in other words, the Gibbs paradox was resolved when we correctly determined the system's initial multiplicity.
We now prove their "equality".
\begin{align*}
\Omega_\text{sys}
    &= \Omega(V_1, N_1)\Omega(V_2, N_2) \times  {N_1+N_2 \choose N_1,N_2} \\
    &= f(N_1)f(N_2)\left(\frac{V_1}{\lambda^3}\right)^{N_1}\left(\frac{V_2}{\lambda^3}\right)^{N_2}\frac{(N_1+N_2)!}{N_1!N_2!}\delta^2
\end{align*}
Take the logarithm, i.e. calculate the entropy.
\begin{align*}
\ln\Omega_\text{sys}
    &= \ln f(N_1) + \ln f(N_2) + N_1\ln\left(\frac{V_1}{\lambda^3}\right) +N_2\ln\left(\frac{V_2}{\lambda^3}\right) + \ln\left(\frac{(N_1+N_2)!}{N_1!N_2!}\right) + 2\ln \delta
\end{align*}
Discard the constant term.
\begin{align*}
    &\approx \ln f(N_1) + \ln f(N_2) + N_1\ln\left(\frac{V_1}{\lambda^3}\right) +N_2\ln\left(\frac{V_2}{\lambda^3}\right) + \ln\left(\frac{(N_1+N_2)!}{N_1!N_2!}\right)
\end{align*}
One can show that $\ln f(N) \sim 3N/2$ as $N\to\infty$.
\begin{align*}
    &\approx \frac{3(N_1+N_2)}{2} + N_1\ln\left(\frac{V_1}{\lambda^3}\right) +N_2\ln\left(\frac{V_2}{\lambda^3}\right) + \ln\left(\frac{(N_1+N_2)!}{N_1!N_2!}\right)
\end{align*}
Recall that $N_1/V_1=N_2/V_2$.
\begin{align*}
    &= \frac{3(N_1+N_2)}{2} + N_1\ln\left(\frac{V_1/N_1}{\lambda^3}\right) +N_2\ln\left(\frac{V_2/N_2}{\lambda^3}\right) + N_1\ln N_1 + N_2\ln N_2 + \ln\left(\frac{(N_1+N_2)!}{N_1!N_2!}\right) \\
    &= \frac{3(N_1+N_2)}{2} + (N_1+N_2)\ln\left(\frac{V_1/N_1}{\lambda^3}\right) + N_1\ln N_1 + N_2\ln N_2 + \ln\left(\frac{(N_1+N_2)!}{N_1!N_2!}\right)
\end{align*}
Use Stirling's approximation.
\begin{align*}
    &= \frac{3(N_1+N_2)}{2} + (N_1+N_2)\ln\left(\frac{V_1/N_1}{\lambda^3}\right) + (N_1+N_2)\ln(N_1+N_2) \\
    &= \frac{3(N_1+N_2)}{2} + (N_1+N_2)\ln\left(\frac{(N_1+N_2)V_1/N_1}{\lambda^3}\right) \\
    &= \frac{3(N_1+N_2)}{2} + (N_1+N_2)\ln\left(\frac{V_1+V_2}{\lambda^3}\right)
\end{align*}
Therefore, for large $N_1$ and $N_2$,
$$
\ln\Omega_\text{sys} \approx \frac{3(N_1+N_2)}{2} + (N_1+N_2)\ln\left(\frac{V_1+V_2}{\lambda^3}\right) \,.
$$
On the other hand, the entropy after the divider is removed is
\begin{align*}
\ln\Omega_\text{sys}'
&= \ln(\Omega(V_1+V_2, N_1+N_2)) \\
&= \ln f(N_1+N_2) + (N_1+N_2)\ln\left(\frac{V_1+V_2}{\lambda^3}\right) + \ln \delta\\
&\approx \frac{3(N_1+N_2)}{2} + (N_1+N_2)\ln\left(\frac{V_1+V_2}{\lambda^3}\right) \\
&= \ln\Omega_\text{sys} \,.
\end{align*}
Therefore, both $\Omega_\text{sys}$ and $\Omega_\text{sys}'$ lead to the same entropy, and there is no entropy gained by removing the divider.
As a side note, we've only shown that $\Omega_\text{sys}$ and $\Omega_\text{sys}'$ are equal when ignoring low order terms. Actually, we should not expect them to be equal exactly either. Before removing the divider, each box contained a definite number of particles, but after the divider is removed, the particle number in each box will be variable.
Thus, the Gibbs paradox is resolved without needing to use correct Boltzmann counting. If the "naive" approach to calculating the multiplicity was correct all along, why should we still insist on dividing by $N!$ to account for indistinguishability? How is division by $N!$ even justified?
There is one downside to not using correct Boltzmann counting, namely, entropy is no longer extensive. By an extensive entropy, I mean, one where the following is true. The entropy of a total system, is equal to the sum of the entropies of 'each subsystem when considered in isolation', regardless if some subsystems are "identical". This is equivalent to requiring that the multiplicity of the total system, be a product of the multiplicities of 'each subsystem when considered in isolation'.
$$
\Omega_\text{total} = \prod_i \Omega_\text{subsystem $i$}
$$
The multiplicity we have been using so far, clearly does not satisfy this property, we had to multiply by a multinomial coefficient to get the correct total system multiplicity;
$$
\Omega_\text{total} \neq \Omega_\text{left box}\times\Omega_\text{right box} \,,
$$
where $\Omega_\text{left box}$ is the multiplicity of the left box, when the left box is considered in isolation, and similarly for $\Omega_\text{right box}$.
There is a natural connection between how we have been calculating multiplicities so far, and correct Boltzmann counting. By rearranging the equation relating the total system multiplicity (before the divider was removed) and the multiplicities of the left and right boxes (when considered in isolation), we find
\begin{gather*}
\Omega_\text{sys}
    = \Omega(V_1, N_1)\Omega(V_2, N_2) \times  \frac{(N_1+N_2)!}{N_1!N_2!} \\
\frac{\Omega_\text{sys}}{(N_1+N_2)!} = \frac{\Omega(V_1, N_1)}{N_1!}\frac{\Omega(V_2, N_2)}{N_2!}
\end{gather*}
Correct Boltzmann counting absorbs the factorials into the multiplicities, so if we defined
\begin{align*}
\Omega_\text{sys}^{(B)} &\equiv \frac{\Omega_\text{sys}}{(N_1+N_2)!} \\
\Omega^{(B)}(V, N) &\equiv \frac{\Omega(V, N)}{N!}
\end{align*}
we restore the familiar
$$
\Omega_\text{sys}^{(B)} = \Omega^{(B)}(V_1, N_1)\Omega^{(B)}(V_2, N_2)
$$
Since $\Omega_\text{sys}^{(B)}$ differs from $\Omega_\text{sys}$ only by a constant multiple (total particle number $N_1+N_2$ is constant here), maximization of entropy defined as $k_B\ln \Omega_\text{sys}^{(B)}$ and maximization of entropy defined as $k_B\ln \Omega_\text{sys}$, will predict the same equilibrium states.
Essentially what we've shown is that the multiplicity of a monatomic ideal gas, calculated with or without correct Boltzmann counting (hereafter CBC), leads to an entropy which when maximized, predicts the same equilibrium states. One might justify the multiplicity calculated without CBC by Liouville's theorem and the ergodic hypothesis (to argue that macrostate probability is proportional to phase space volume). CBC multiplicity can then be rationalized as being empirically equivalent to the non-CBC multiplicity. Similar arguments could probably be made for legitimizing the use of CBC in other circumstances where CBC might not be "easily justified" from the foundational physical assumptions.

Why do we care about an extensive entropy?
I think the reason for this is that classical thermodynamics assumes that entropy is an extensive quantity, i.e. the entropy of the total system can always be calculated as the sum of entropies of each subsystem when considered in isolation.
For example, consider the system above in the setup for Gibbs paradox, and suppose $N_1/V_1\neq N_2/V_2$. We wish the find the equilibrium state when the divider is removed. Let $S_\text{sys}$ be the total system's entropy, and $S_1$ and $S_2$ be entropies of the left and right box subsystems (when considered in isolation) respectively. Assume entropy is extensive.
\begin{align*}
\frac{\partial S_\text{sys}}{\partial N_1}
&= \frac{\partial S_1}{\partial N_1} + \frac{\partial S_2}{\partial N_1}\\
&= \frac{\partial S_1}{\partial N_1} - \frac{\partial S_2}{\partial N_2} \\
\end{align*}
At equilibrium, $\partial S_\text{sys}/\partial N_1=0$.
\begin{align*}
\frac{\partial S_1}{\partial N_1} &= \frac{\partial S_2}{\partial N_2} \\
\mu_1 &= \mu_2
\end{align*}
If entropy is not extensive, this will not be true. Without correct Boltzmann counting, the total system multiplicity is (where we use $\Omega(E, V, N)$ so partial derivatives with respect to $N$ are when $E$ and $V$ are held constant, and not $T$ and $V$)
\begin{gather*}
\Omega_\text{sys}
= \Omega(E_1, V_1, N_1)\Omega(E_2, V_2, N_2)\times\frac{(N_1+N_2)!}{N_1!N_2!} 
\end{gather*}
The entropy is then, where $S_i = k_B\ln \Omega(E_i, V_i, N_i)$,
\begin{gather*}
S_\text{sys} = S_1 + S_2 - k_BN_1\ln N_1 - k_BN_2\ln N_2 + k_B\ln((N_1+N_2)!) - k_B(N_1+N_2)
\end{gather*}
Differentiate with respect to $N_1$, keeping $N_1+N_2$ constant.
\begin{gather*}
\frac{\partial S_\text{sys}}{\partial N_1} = \frac{\partial S_1}{\partial N_1} - \frac{\partial S_2}{\partial N_2} - k_B\ln N_1 + k_B\ln N_2
\end{gather*}
At equilibrium, $\partial S_\text{sys}/\partial N_1 = 0$.
\begin{gather*}
\frac{\partial S_1}{\partial N_1} - k_B\ln N_1 = \frac{\partial S_2}{\partial N_2} - k_B\ln N_2 \\
\mu_1 + k_BT\ln N_1 = \mu_2 + k_BT\ln N_2
\end{gather*}
This is a different (equally valid) equilibrium condition than was found before.
I don't think it is possible to "test" if gas particles in the real world have extensive or intensive entropies by "measuring their chemical potentials", since (my understanding is) when we "measure the chemical potential" of some chemical species, we are really only fitting parameters to a thermodynamic model that already assumes an extensive entropy, so using this method to prove that "entropy of real gases are extensive" might be circular and meaningless.
How do we justify correct Boltzmann counting for chemical interactions between subsystems of different chemical species, or when particle numbers are not conserved?
This is a question that came to my mind that others might have to. Ultimately, I think this question is kind of meaningless. (I think) the thermodynamic potentials of real gases are determined through experiments, and "parameter fitting", and not through "first principles" as we have done for the ideal monatomic gas.
The chemical potential determined by the Sackur-Tetrode equation might not be comparable to chemical potentials of real gases determined through experiments, since these quantities are (probably?) relative to some "standard zero" potential, whereas it's not clear what the chemical potential determined through the Sackur-Tetrode equation is relative to.
Also, Liouville's theorem is no longer applicable if total particle numbers (of some particles that evolve via Hamilton's equations) is not conserved, so correct Boltzmann counting is no less, or more, easier to justify than "incorrect Boltzmann counting". We would need another way to reason as to why the microstates we have chosen ought to be equiprobable.
As a side note, the goal of statistical mechanics (seems to be) to use probabilistic arguments to determine thermodynamic potentials of new systems, such that they interact "consistently" with the existing thermodynamic potentials of studied systems.
A cool observation.
As another side note, going back to the Gibbs paradox example, when we calculated the multiplicity, without correct Boltzmann counting, after the divider is removed, we could have done it in two ways, either
$$
\Omega_\text{sys}' = \Omega(V_1+V_2, N_1 + N_2) \,.
$$
as was done initially, or
$$
\Omega_\text{sys}' = \Omega(V_1+V_2, N_1)\Omega(V_1+V_2, N_2) \,.
$$
The two can be shown to be equivalent (again, not exactly equivalent since the second calculation does not consider the microstates where a small group of particles, a single particle say, contains all the energy of the system, and the others are very slow moving). The first interpretation says, the system after the divider is removed, is the same as a box of volume $V_1+V_2$ with $N_1+N_2$ particles. The second interpretation says, the system after the divider is removed, is the same as expanding the volume of the left and right boxes, to fill the size of both boxes, and the left and right boxes kind of "coexist".
The second interpretation is also how one would calculate the total system multiplicity for distinct gases that mix after the divider is removed. Without correct Boltzmann counting, the multiplicities of two mixed distinct gases, is the same as two mixed identical gases, and the multiplicity of two unmixed distinct gases, is lower than the multiplicity of two unmixed identical gases. This is opposite to correct Boltzmann counting. In some sense, this is more intuitive. Two mixed distinct gases is just as "useless" as two mixed identical gases, but useful work can be extracted from two unmixed distinct gases, so perhaps its entropy "should" be lower than two "unmixed/mixed" identical gases.

Sorry if this was too long, I just thought some parts might be interesting to others. I'm still learning the subject as well so I don't know if everything said is strictly correct, but I'm open to feedback.
This paper
Jaynes, E. T. (1996). "The Gibbs paradox".
by Jaynes was also kind of interesting as well. In one part he talks about how entropy is only defined up to an abitrary function of particle number $N$ in "phenomenological thermodynamics".

A:

There is no need for quantum mechanics. You can find a classical argument in Landau-Lifshitz, Statistical Physics, Chapter III, Paragraph 31.
What you have to understand is that when integrating over phase space, for example to calculate the free energy in the canonical ensemble, you should not integrate over all values of the canonical variables $p,q$, but only over the values which correspond to different physical states. Landau uses the notation $\int'$ to denote this; for example for the free energy you will have
$$\beta F = -  \ln \int' e^{-\beta \mathcal H(p,q)} dp dq$$ 
where $\beta=1/k_B T$ as usual. 
Now, it is often more practical to just integrate over the whole phase, i.e. to take $\int dp dq$ space instead that over the portions of phase space corresponding to different states, i.e. $\int' dp dq$. 
However, if you do so, you are overcounting the number of states. In the case of a gas of $N$ atoms, in order to get the correct counting, you have to divide the integral extended over the whole phase space by the number of possible permutations of $N$ atoms, which is $N!$. This is because if you switch the position of two identical atoms you are not getting a new physical state (this is the key point of the whole argument). To sum up, for a gas of $N$ identical atoms we have
$$\int ' \dots dp dq = \frac 1 {N!} \int \dots dp dq$$
which is Eq. 31.7 in my edition of Landau-Lifshits.

