Proof of Coulomb's law in two and higher dimensions I have found that the Coulomb force in two dimension varies with $\frac 1 r$:
\begin{equation}\tag{2}F=\frac{1}{2\pi\epsilon}\cdot\frac{q_1q_2}{r}\end{equation}
But I was not able to prove it. I think it can be proved using the Laplace equation and other further modifications.
One question related to this topic is answered here:

2 dimensional Coulomb's law equation

But the proof has not been given, rather a general introduction is given. Also, some links related to Green's function have been provided, which I didn't understand. (If anyone can elaborate on this, that would be helpful)
so my question is
How to prove Coulomb's law in two dimensions or genralised N dimension?
And  Coulomb's law depends upon $ \vec{r} $ only and,
\begin{equation}\tag{2}F=\frac{1}{4\pi\epsilon}\cdot\frac{q_1q_2}{r^2}\end{equation}
so we can put the value of $ r^2 $ (distance) both in 2 or 3 dimensions, so why would the equation change as it is not dependent upon $\Theta $ (Theta) or $\Phi $ (Phi). So intuitively the equation should not change?
Is Gauss law applicable only in 3 dimensions or valid for any dimension, because the Gauss divergence theorem is only for 3 dimensions?
Also, please write some other insightful details if comes across.
 A: The strict answer to "How does electrostatics work in higher dimensions?" is "Nobody knows", because we can't exactly pop over to a 5-D universe and do experiments.  So if you want to theorize about how physical laws would work in higher dimensions, you basically have to write down our Universe's laws in a way that generalizes in a simple way to these higher dimensions.
For electrostatics, Gauss's Law (in differential form) generalizes in a simple way:  we can write down the divergence of a higher-dimensional vector field as
$$
\vec{\nabla} \cdot \vec{E} = \frac{\partial E_1}{\partial x_1} + \frac{\partial E_2}{\partial x_2} + \frac{\partial E_3}{\partial x_3} + \dots = \frac{\rho}{\epsilon_0},
$$
where $\rho$ is now charge per unit volume in $N$ dimensions.
This can then be shown mathematically to be equivalent to saying that
$$
\oint \vec{E} \cdot d^{N-1} \vec{a} = \frac{Q_\text{enc}}{\epsilon_0},
$$
where the integral on the left-hand side is over an $N-1$ dimensional surface, and $Q_\text{enc}$ is the amount of charge enclosed in that surface.
We can similarly define a notion of "spherical symmetry" in these higher dimensions. Assuming that the field of a point charge is spherically symmetric in these higher dimensions, we can choose a "sphere" of radius $r$ to integrate over for Gauss's Law, with the result that
$$
|\vec{E}(r)| A_N(r) = \frac{Q_\text{enc}}{\epsilon_0},
$$
where $A_{N-1}(r)$ is the "surface area" of an $N-1$-dimensional "sphere" of radius $r$.  These surface areas can be calculated, with the result that
$$
A_{N-1}(r) = \frac{2 \pi^{N/2}}{\Gamma\left(\frac{N}{2}\right)} r^{N-1}.
$$
Thus, given the above assumptions, the field of a point charge in $N$ dimensions should be
$$
|\vec{E}(r)| = \frac{\Gamma\left(\frac{N}{2}\right)} {2 \pi^{N/2} \epsilon_0}\frac{q}{r^{N-1}}.
$$
A: As with all derivations, it depends on what you want to treat as fundamental. Typically we would derive Coulomb's law from the Maxwell equations, so we're trying to solve
$$\nabla\cdot \mathbf{E} = -\nabla^2 \varphi = q\delta(\mathbf{x})/\epsilon_0\qquad (1)$$
In $n$ spatial dimensions and in Cartesian coordinates $(x_1,\ldots,x^n)$, this becomes
$$\sum_{k=1}^n \frac{\partial^2}{\partial x_k^2} \varphi = -\frac{q}{\epsilon_0}\delta(\mathbf x)\qquad(2)$$
Because this problem has spherical symmetry, we can move to hyperspherical coordinates.  If we do so, we will find that$^\dagger$
$$\frac{1}{r^{n-1}}\frac{\partial }{\partial r}\left(r^{n-1} \frac{\partial\varphi}{\partial r}\right) = -\frac{q}{\epsilon_0} \delta(r)\qquad (3)$$
Away from $r=0$, we would therefore have that $$\frac{\partial}{\partial r}\left(r^{n-1} \frac{\partial \varphi}{\partial r}\right)=0 \implies r^{n-1} \frac{\partial \varphi}{\partial r} = c$$ for some constant $c$, and therefore that $\varphi = c\  r^{2-n}+d$ (unless $n=2$, in which case we'd have a logarithm).  The constant $d$ can be set to zero by demanding that the potential vanish at infinity (this is an arbitrary choice, but a convenient one).  The constant $c$ can be determined by using the divergence theorem to integrate $(1)$ over a hypersphere of radius $R$.  Because of the spherical symmetry, the left hand side would be the surface area of the $(n-1)$-sphere of radius $R$ times $\varphi'(R)$:
$$-\frac{2\pi^{n/2}}{\Gamma(n/2)}R^{n-1} \varphi'(R)=\left(\frac{2\pi^{n/2}(n-2)}{\Gamma(n/2)}\right) c$$
while the right hand side is simply equal to $q/\epsilon_0$ because of the delta function.  As a result,
$$\varphi(r) = \frac{\Gamma(n/2)}{2(n-2)\pi^{n/2}\epsilon_0} \frac{q}{r^{n-2}}\qquad (4)$$
In $n=3$ dimensions, we have $\Gamma(3/2)=\sqrt{\pi}/2$ so this reduces to the familiar case
$$\varphi^{(3)}(r) = \frac{1}{4\pi\epsilon_0} \frac{q}{r} \implies \mathbf{E}^{(3)}(r) = \frac{1}{4\pi\epsilon_0} \frac{q}{r^2}\hat r$$
In 4-dimensions, $\Gamma(2)=1$ so we would have
$$\varphi^{(4)}(r) = \frac{1}{4\pi^2 \epsilon_0} \frac{q}{r^2} \implies \mathbf{E}^{(4)}(r) = \frac{1}{2\pi^2 \epsilon_0} \frac{q}{r^3} \hat r$$
In the other direction, for $n=1$ we have $\Gamma(1/2)=\sqrt{\pi}$ and so
$$\varphi^{(1)}(r) = -\frac{1}{2\epsilon_0} q r \implies \underbrace{\mathbf{E}^{(1)}(r)=\frac{1}{2\epsilon_0} q \hat r}_{\text{constant}}$$


So intuitively the equation should not change?

The problem is that $\nabla^2$ changes in higher dimensions, so if you reuse the familiar form of Coulomb's law then it will not obey the Maxwell equations. Assuming you'd like to treat the latter as more fundamental, we need to use Gauss' law to find the more general form of Coulomb's law.

Is Gauss law applicable only in 3 dimensions or valid for any dimension, because the Gauss divergence theorem is only for 3 dimensions?

The divergence theorem holds in an arbitrary number of dimensions.  If we assume that Gauss' law holds in an arbitrary number of dimensions, then we find Coulomb's law as I did above.  Of course, Gauss' law is a physical statement, not a purely mathematical one, so there's no way to mathematically prove that it holds for all dimensions.

$^\dagger$This expression should not be taken too literally, as the delta function at the origin has some pathological issues in spherical coordinates.  The spirit of this equation is that we will find the solution for $r\neq 0$, and obtain the remaining undetermined constant by integrating $(1)$.
A: In addition to the great answer above.
Explaining intuitively how the Coulomb law looks in n dimensions.
In two dimensions:
$$F=\frac 1 {2\pi\epsilon}\frac  q r$$
In three dimensions:
$$F=\frac 1 {4\pi\epsilon}\frac{q}{r^2}$$
You are familiar with this. In two dimensions, the force varies linearly with $r$. So the force is "diluted" and varies linearly with $\frac 1 r$.
The same reasoning holds for every higher dimension higher. The force field lines get distributed over a spherical hyper-volume.
So:
$$F_n=\frac{q}{V_{n-1}\epsilon},$$
where $V_{n-1}$ is the area of an $(n-1)$-sphere (the derivatine wrt $r$ of the formula of the volume of an $n$-sphere), to be found in this article.

Does spherical symmetric assumption would be valid in N dimensions or is it controversial?

Due to string theoretical considerations, it could be that on a small scale the inverse square law doesn't hold anymore.
Just as is the case (in the string theoretical approach) with the inverse square law of gravity. On a small scale, again due to a small extra space dimension, it could be that the inverse square law changes in an inverse $r^3$ law, though I'm not sure if the symmetry is exactly spherical. If the extra space dimension resembles a rolled-up cylinder (which, according to the modern theory can have observable consequences because the radius of the cylinder is well above the Planck length).
