Short Answer: Every particle attracts every other, so the total attraction (which is the correction to pressure) has to be proportional to $N^2$.

Long Answer: If you want to rigorously derive the VdW equation of state, you should look into the Cluster expansion. It is covered in most textbooks on statistical mechanics. For a good introduction, have a look at David Tong's lecture notes.

More intuitively, you can view the VdW equation as a mean field theory of a weakly interacting gas. Take a look at the partition function $$ Z = \frac{1}{N! \lambda^{3N}} \int e^{-\beta \sum U(\Delta x_{ij})} dx^{3N} \,, $$ with the momentum integration already performed and put into the thermal wavelength $\lambda$.

For an arbitrary potential $U$, this integral is impossible to do exactly. But you can do a mean field approximation, where you assume that all particles feel the same 'mean field potential' $U_{\rm MF}(r)$. The partition function then factorizes into $$ Z = \frac{1}{N! \lambda^{3N}} \left[\int e^{-\beta U_{\rm MF}(x)} dx^3\right]^N \,. $$

What can you say about $U_{\rm MF}(r)$? You know two things:

1. Particles can't overlap, so that $U_{\rm MF}(r < R) = \infty$. You therefore have an exlcluded volume of size $nb$, where $b$ is you particle volume. This gives the correction to the $V$ term in the VdW equation.
2. The potential will be slightly attractive at long range. Since every particle attracts every other, $U_{\rm MF}(r >> R)$ should be proportional to the particle density, say $U_{\rm MF} = a N / V$.

Since our simplified $U_{\rm MF}$ does not depend on position anymore, we can pull it out of the integral. $$ Z = \frac{(e^{\beta a N / V})^N}{N! \lambda^{3N}} \left[\int dx^3\right]^N = \frac{1}{N!} \left[\frac{e^{\beta a N / V} (V- Nb)}{\lambda^3}\right]^N\,. $$ Here you can see that the $N^2$ dependence comes from the overall exponent of $N$: $$ \left(e^{\beta a N / V}\right)^N = e^{\beta a N^2 / V} $$ You can derive the VdW equation of state by calculating the Gibbs free energy $G$ and looking at its derivatives. Along the way you will pick up another $V$ to give you $a N^2 / V^2$.

If you want to know more I recommend the book Introduction to the Theory of Soft Matter by Jonathan Selinger.

Short Answer: Every particle attracts every other, so the total attraction (which is the correction to pressure) has to be proportional to $N^2$.

Long Answer: If you want to rigorously derive the VdW equation of state, you should look into the Cluster expansion. It is covered in most textbooks on statistical mechanics. For a good introduction, have a look at David Tong's lecture notes.

More intuitively, you can view the VdW equation as a mean field theory of a weakly interacting gas. Take a look at the partition function $$ Z = \frac{1}{N! \lambda^{3N}} \int e^{-\beta \sum U(\Delta x_{ij})} dx^{3N} \,, $$ with the momentum integration already performed and put into the thermal wavelength $\lambda$.

For an arbitrary potential $U$, this integral is impossible to do exactly. But you can do a mean field approximation, where you assume that all particles feel the same 'mean field potential' $U_{\rm MF}(r)$. The partition function then factorizes into $$ Z = \frac{1}{N! \lambda^{3N}} \left[\int e^{-\beta U_{\rm MF}(x)} dx^3\right]^N \,. $$

What can you say about $U_{\rm MF}(r)$? You know two things:

1. Particles can't overlap, so that $U_{\rm MF}(r < R) = \infty$. You therefore have an exlcluded volume of size $nb$, where $b$ is you particle volume. This gives the correction to the $V$ term in the VdW equation.
2. The potential will be slightly attractive at long range. Since every particle attracts every other, $U_{\rm MF}(r >> R)$ should be proportional to the particle density, say $U_{\rm MF} = -a N / V$.

Since our simplified $U_{\rm MF}$ does not depend on position anymore, we can pull it out of the integral. $$ Z = \frac{(e^{\beta a N / V})^N}{N! \lambda^{3N}} \left[\int dx^3\right]^N = \frac{1}{N!} \left[\frac{e^{\beta a N / V} (V- Nb)}{\lambda^3}\right]^N\,. $$ Here you can see that the $N^2$ dependence comes from the overall exponent of $N$: $$ \left(e^{\beta a N / V}\right)^N = e^{\beta a N^2 / V} $$ You can derive the VdW equation of state by calculating the Gibbs free energy $G$ and looking at its derivatives. Along the way you will pick up another $V$ to give you $a N^2 / V^2$.

If you want to know more I recommend the book Introduction to the Theory of Soft Matter by Jonathan Selinger.