Screening of a coulomb-like potential

Imagine that we have a potential in the form of: $$U(r) \propto \frac{1}{r^n}$$ in a 2D system, with the high concentration of particles interacting with the above potential. How do you find the screening length of this system? I am asking this question because I need to consider a cut-off radius for a simulation and considering such a cut-off radius could enhance my simulation speed a lot. I want to make sure that I am choosing the right value for that.

Thanks

Linearizing a Poisson-Boltzmann type equation (typically expressed in terms of different types of charges in solution, here only using one such term; $$e$$: elementary charge, $$c_0$$: concentration constant, $$\epsilon$$: relative permittivity, $$\epsilon_0$$: permittivity of vacuum, $$k_{\rm B}$$: Boltzmann constant, and $$T$$: temperature): $$\nabla^2 \phi = \frac{c_0\,e}{\epsilon\,\epsilon_0}\exp\left(-\frac{e\phi}{k_{\rm B}T}\right) \approx \frac{c_0\,e}{\epsilon\,\epsilon_0} \left( 1 - \frac{e\phi}{k_{\rm B}T} \right),$$ the screening term can be identified as the additional term applied to the potential $$\phi$$ besides $$\nabla^2$$: $$\left( \nabla^2 + \frac{c_0\,e^2}{\epsilon\,\epsilon_0\,k_{\rm B}T} \right) \phi = \frac{c_0\,e}{\epsilon\,\epsilon_0}.$$
Fourier transformation of $$U_{\rm screened}(r) \propto r^{-n} \exp(-k\,r)$$ will allow for determining $$k$$ by comparison to the above equation. In 3D with the usual Coulomb interaction $$\propto\!r^{-1}$$, e.g., one would have for the Fourier transform: $$\propto\!q^{-2} + k^{-2}$$ (where $$q$$ is the wave vector and $$k$$ is the inverse screening length). $$q^{-2}$$ can be identified with $$\nabla^2$$ in Fourier space, and $$k^{-1}_{{\rm 3D},n=1} \propto \sqrt{\frac{c_0\,e^2}{\epsilon\,\epsilon_0\,k_{\rm B}T}} .$$
For 2D, Fourier transformation of $$U_{\rm screened}(r)=r^{-n}\exp(-kr)$$ yields the following: $$\int_0^{2\pi} {\rm d}\psi \int_0^\infty {\rm d}r \exp[iqr\cos(\psi)] \, r \, U_{\rm screened}(r) = 2\pi \int_0^\infty J_0(qr)\,r\,U_{\rm screened}(r) \, {\rm d}r = \frac{2\pi}{\sqrt{q^2+k^2}} ,$$ for $$n\ge 1$$. $$J_0$$ is the Bessel function of first kind and order zero. Due to the square root, unfortunately no direct comparison to the linearized Poisson-Boltzmann equation can be made. For $$q\rightarrow 0$$ (corresponding to long wavelengths in real space), $$\frac{1}{\sqrt{q^2+k^2}}\rightarrow k^{-1}$$ and thus $$k^{-1}_{\rm 2D,longrange} \approx \frac{c_0\,e^2}{2\pi\,\epsilon\,\epsilon_0\,k_{\rm B}T} ,$$ interestingly not dependending on the value of $$n\ge 1$$.