  In the past few decades, much attention has been paid to the bubbling problem
for semilinear Neumann elliptic equation with the critical and subcritical
polynomial nonlinearity, much less is known if the polynomial nonlinearity is
replaced by the exponential nonlinearity. In this paper, we consider the
following semilinear Neumann elliptic problem with the Trudinger-Moser
exponential growth: \begin{equation*}\begin{cases} -d\Delta
u_d+u_d=u_d(e^{u^2_{d}}-1)\ \ \mbox{in}\ \Omega,\\ \frac{\partial
u_d}{\partial\nu}=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{on}\
\partial \Omega,\\ \end{cases}\end{equation*} where $d>0$ is a parameter,
$\Omega$ is a smooth bounded domain in $\mathbb{R}^2$, $\nu$ is the unit outer
normal to $\partial \Omega$. We first prove the existence of a ground state
solution to the above equation. If $d$ is sufficiently small, we prove that any
ground state solution $u_d$ has at most one maximum point which is located on
the boundary of $\Omega$ and characterize the shape of ground state solution
$u_d$ around the condensation point $P_d$. The key point of the proof lies in
proving that the maximum point $P_d$ is close to the boundary at the speed of
$\sqrt{d}$ when $d\rightarrow0$ and $u_d$ under suitable scaling transform
converges strongly to the ground state solution of the limit equation $\Delta
w+w=w(e^{w^2}-1)$.
  Our proof is based on the energy threshold of cut-off function, the
concentration compactness principle for the
  Trudinger-Moser inequality, regularity theory for elliptic equation and an
accurate analysis for the energy of the ground state solution $u_d$ as
$d\rightarrow0$. Furthermore, by assuming that $\Omega$ is a unit disk, we
remove the smallness assumption on $d$ and show the maximum point of ground
state solution $u_d$ must lie on the boundary of $\Omega$ for any $d>0$.
