Let $E_{\tau}:= \mathbb{C}/(\mathbb{Z}+ \mathbb{Z} \tau)$ be a flat torus and $G(z; \tau)$ be the Green function on $E_{\tau}$. Consider the multiple Green function $G_{n}$ on$(E_{\tau})^{n}$: \[ G_n (z_1, \cdots ,z_n ; \tau) := \sum_{i \lt j} G(z_i - z_j ; \tau) - n \sum_{i=1}^n G(z_i ; \tau). \] We prove that for $ \tau \in i \mathbb{R}_{\gt 0}$, i.e. $E_\tau$ is a rectangular torus, $G_n$ has exactly $2n + 1$ critical points modulo the permutation group $S_n$ and all critical points are non-degenerate. More precisely, there are exactly $n$ (resp. $n+1$) critical points $ \boldsymbol{a}$’s with the Hessian satisfying $(-1)^n \det D^2 G_n (\boldsymbol{a} ; \tau) \lt 0$ (resp. $\gt 0$). This confirms a conjecture in [4]. Our proof is based on the connection between $G_n$ and the classical Lamé equation from [4, 19], and one key step is to establish a precise formula of the Hessian of critical points of $G_{n}$ in terms of the monodromy data of the Lamé equation. As an application, we show that the mean field equation \[ \Delta_u + e^u = \rho \delta_0 \textrm{ on } E_\tau \] has exactly $n$ solutions for $8 \pi n - \rho \gt 0$ small, and exactly $n+1$ solutions for $\rho - 8 \pi n \gt 0$ small.