Why is the density of the Fermi gas in a neutron star not changing the potential depth caused by the strong nuclear interaction?

You are quite correct that a neutron star is not supported by ideal neutron degeneracy pressure. Any book or web source that claims so should be given a wide berth.

As long ago as 1939 Oppenheimer & Volkhoff showed that a neutron star supported by ideal NDP became unstable at finite density, with a maximum mass of around $0.7 M_{\odot}$. All measured neutron star masses are much higher than this.

The repulsive core of the strong nuclear force in asymmetric nuclear matter is almost certainly what supports neutron stars. The polytropic index of the pressure can exceed 2, as opposed to somewhere between 4/3 and 5/3 for ideal NDP, so is a much harder equation of state.

The review by Lattimer (2013) does a good job of describing how observations of neutron star masses and radii provide constraints on the equation of state and uncertain parameters in the symmetry energy beyond nuclear densities.

My answer concerns the origin of the "repulsive hard core" referred to in @RobJeffries answer (with which I totally concur). This mechanism is responsible for the relative constancy of the central nuclear density of large nuclei as well as the stability of neutron stars below a certain mass limit. If you have read my answer to this question: Why is nuclear force spin dependent? you will see that within the relativistic mean field (or Hartree) approach to the nuclear many-body problem the central nuclear potential has two terms: 1) $U_s$, an attractive term arising from the exchange of the sigma (scalar, isoscalar) meson, and 2) $U_0$, a repulsive term arising from the exchange of the omega (vector, isoscalar) meson. These empirically determined potentials are the result of an elaborate self-consistent field (Hartree or Hartree-Fock) procedure that involves solving by numerical integration the Dirac equation $$[c\vec{ \alpha } \cdot \vec{ p}+\gamma^0(m+U_s+\gamma^0U_0)]\psi=E\psi$$for single particle orbitals for doubly magic nuclei over the entire periodic table. The wave function $\psi$ is a Dirac 4-component spinor that takes the form: $$\psi(\vec r)=\frac{1}{r}\binom {F(r)Y^{\omega}_{jm}(\hat r)}{iG(r)Y^{-\omega}_{jm}(\hat r)}$$ when the potentials are spherically symmetric. The functions $F(r)$ and $G(r)$ are the large and small component radial wave functions and $Y^{\omega}_{jm}(\hat r)$ is the Pauli central field spinor. Once a set of trial orbitals for a particular nucleus are found, the nuclear potentials are obtained by folding the relevant nucleon density functions with the associated Yukawa function. Iteration continues until two successive sets of orbitals and potentials are found to be in agreement (self-consistency).

The masses and couplings of these two mesons are the primary parameters in the model. When the work was done initially only the omega mass was known so the other parameters were adjusted to yield agreement with the experimental binding energies, charge distributions, and single particle separation energies of these doubly magic nuclei. The single particle eigenvalues turned out to have the correct ordering needed to reproduce the known magic number sequences for neutron and proton shells. This was a result of the large spin-orbit interaction inherent in this combination of attractive and repulsive potentials. The nuclear potentials resulting from this empirical procedure can be roughly described as follows: the attractive (scalar) potential ($U_s(r)$) has a depth between 500-600 MeV while the repulsive (vector) potential ($U_0(r)$) has a height between 400-500 MeV. In a zeroth order non-relativistic approximation the nuclear potential would be the sum of these two ($U_s(r)+U_0(r)$, a well with a depth of 50-100 MeV).

So what negates the argument:

"At least the (attractive) Hartree-Potential in the Hartree-Fock Formalism grows linearly with nulceon density (which grows by $\frac{1}{V}$). This gain in $\frac{1}{V}$ would win over the (repulsive) increase in the kinetic energy under compression, which grows only by $V^{-\frac{2}{3}}$). Thus, the negative pressure caused by the increasing binding energy grows faster than the positive Fermi gas pressure due to the kinetic energy U. It is energetically advantageous to collapse."

and prevents collapse (even in finite nuclei as well as neutron stars)? The answer comes from studying the nucleon densities that are folded with Yukawa functions to yield the scalar ($U_s$) and vector ($U_0$) potentials. The vector nucleon density is the same as the probability density so the above argument applies, except that the vector potential is repusive, not attractive. The scalar density, on the other hand, differs from the probability density (the square of the small component of the radial wave function is subtracted from the square of the large component ($F^2-G^2$) rather than added ($F^2+G^2$) in the case of the vector density. This means that if the nucleon density is squeezed beyond the normal equilibrium value, the repulsive potential indeed grows but the attractive scalar potential will not grow as much since the squeezing will also enhance the size of the small component ($G$) relative to the large component ($F$) ($G$ is the order of $\frac{v}{c}$ compared to F). This is the heart of the nuclear saturation mechanism in this approach. It is more a subtle relativistic effect that enhances the repulsive omega meson exchange relative to the attractive sigma meson exchange as the nucleon density increases. In the end a very stiff EOS results.

References: This work was my dissertation research. The first published calculation (Hartree approximation) was Phys. Rev C5 (1972) 241. The full Hartree-Fock calculation was Phys Rev C9 (1974) 537. My calculations were restricted to finite nuclei. Walecka published results with a similar model for infinite nuclear matter in the Fermi gas approximation (for balanced as well as unbalanced nuclear matter, applicable to neutron stars): Ann of Physics (NY) 83 (1974) 491.