\section{Introduction}

Let $G$ be a finite group, $w(x_1,...,x_n)$ a group word and $w$ the associated word map $w:G^{(n)}\rightarrow G$. For $g\in G$ denote by $N_w^G(g)$ the number of solutions to $w= g$ and by $P_w^G(g)$ the probability that a random $n$-tuple $\textbf{g}=(g_{1},...,g_{n})\in G^{(n)}$ satisfies $w(\textbf{g})=g$, i.e.\
$$
P^G_w(g) = \frac{N^G_w(g)}{|G|^n}.
$$
The study of $P_w^G$ has attracted a lot of attention in recent years. For example, in \cite{DPSSh}, it was shown that for $1\neq w$ and a finite simple group $G$, $P_w^G(1)\rightarrow 1$ as $|G|\rightarrow\infty$. Results in \cite{LaSh} provide sharp bounds on $P_w^G$ for general words $w$ and finite simple groups $G$. If $G$ is abelian, then the word map
$$
w:G^{(n)}\rightarrow G,
$$
is a homomorphism and it is clear that
$$
N_w^G(1)=|\mbox{Ker } w|=\frac{|G|^{n}}{|\mbox{Im }w|}\geq |G|^{n-1}
$$
and so $P_w^G(1)\geq\frac{1}{|G|}$. It is a conjecture of Alon Amit (see \cite{Abert}) that if $G$ is a nilpotent group then $P_w^G(1)\geq\frac{1}{|G|}$. In \cite{Levy1} we prove Amit's Conjecture in the special case where the nilpotency class is $2$. Note that since the statistics $N_w^G(1)$ and $P_w^G(1)$ are multiplicative under direct products, we may reduce to the case where $G$ is a finite $p$-group. In this paper we study the fibre sizes of a particular class of words over $p$-groups. In Section \ref{class2} we show that this is enough to determine the fibre sizes of all word maps over $p$-groups of class $2$ with exponent $p$.

For $t\in\mathbb{N}$, let $c_t$ denote the word map given by $c_t = [x_1,y_1]...[x_t,y_t]$ and, for $g\in G$, let $N_t^G(g)$ denote $N_{c_t}^G(g)$. Similarly let $P_t^G$ denote the \textit{$c_t$-distribution} on $G$, i.e.\
$$
P_t^G(g) = \frac{N_t^G(g)}{|G|^{2t}},
$$
and let $U^G$ be the uniform distribution on $G$ (i.e. $U^G(g) = 1/|G|$). By a classical result of
Frobenius from 1896 (see, for example, \cite{Isaacs}) we have
\begin{equation}\label{equation1}
N_t^G(g) = |G|^{2t-1}\sum_{\chi\in\Irr(G)}\frac{\chi(g)}{\chi(1)^t},
\end{equation}
where $1$ is the identity element of $G$ and we sum over the irreducible complex characters of $G$. In the main result of this paper, Theorem \ref{main}, we give explicit formulae to compute the numbers $N^G_t(g)$, for finite $p$-groups $G$ of nilpotency class less than $p$ with exponent $p$, in terms of the number of rational points of certain algebraic varieties.

In \cite{GarionShalev} Garion \& Shalev prove the following.

\begin{proposition}[Proposition 1.1 \cite{GarionShalev}]
Let $G$ be a finite group. Then
$$
||P_1^G-U^G||_1 \leq  \left(\sum_{\chi\in\Irr(G), \chi\neq\text{Id}}\chi(1)^{-2}\right)^{1/2}
$$
where $||P_1^G-U^G||_1=\sum_{g\in G}|P_1^G(g)-U^G(g)|$
\end{proposition}

As the authors remark, this bound has no content when the sum on the right hand side is greater than or equal to $1$. Since the non-trivial linear characters of $G$ contribute $|G/G'|^{-1}$ to the sum the result can only be useful for perfect groups. Since the maps $c_t$ take values in $G'$ it is not hard to adapt their proof and deduce the following.

\begin{proposition}\label{propdist}
Let $G$ be a finite group. Then
$$
||P_t^{G'}-U^{G'}||_1 \leq \left(\frac{|G'|}{|G|}\sum_{\chi\in\Irr(G), \chi(1)\neq 1}\chi(1)^{-2t}\right)^{1/2}
$$
\end{proposition}

If the right hand side in Proposition \ref{propdist} is close to zero, then the fibres of the maps $c_t$ are roughly the same size and the values are uniformly distributed. The proof follows easily from the following lemmas.
\begin{lemma}\label{prevlemma}
Let $G$ be a finite group. Then
$$
\sum_{g\in G'} P^G_t(g)^2= \frac{1}{|G'|}+\frac{1}{|G|}\sum_{\chi\in\Irr(G),\chi(1)\neq 1}\chi(1)^{-2t}.
$$
\end{lemma}
\begin{proof}
The proof follows from Lemma 2.1 in \cite{GarionShalev} and noting that $P^G_t(g) = 0$ for $g\not\in G'$.
\end{proof}
\begin{lemma}
Let $G$ be a finite group. Then
$$
\sum_{g\in G'}\left(P^G_t(g)-\frac{1}{|G'|}\right)^2=\frac{1}{|G|}\sum_{\chi\in\Irr(G),\chi(1)\neq 1}\chi(1)^{-2t}.
$$
\end{lemma}
\begin{proof}
By Lemma \ref{prevlemma},
\begin{eqnarray*}
\sum_{g\in G'}\left(P^G_t(g)-\frac{1}{|G'|}\right)^2 &=& \sum_{g\in G'}P^G_t(g)^2-\frac{2}{|G'|}\sum_{g\in G'}P^G_t(g)+\frac{1}{|G'|}\\ &=&\frac{1}{|G|}\sum_{\chi\in\Irr(G),\chi(1)\neq 1}\chi(1)^{-2t}
\end{eqnarray*}
since $\sum_{g\in G'}P^G_t(g) = 1$.
\end{proof}

\begin{proof}[Proof of Proposition \ref{propdist}]
This follows from the Cauchy-Schwarz inequality,
$$
(||P_t^{G'} - U^{G'}||_1)^2 = \left(\sum_{g\in G'}(P_t(g)-\frac{1}{|G'|}\right)^2\leq|G'|\sum_{g\in G'}\left(P_t(g)-\frac{1}{|G'|}\right)^2
$$
and the previous lemma.
\end{proof}

In Section \ref{mainsec} we will develop formulae for the fibre sizes of the word maps $c_t$ over $p$-groups of nilpotency class less than $p$ with exponent $p$. Moreover, these results extend, more generally, to $p$-groups obtained by `base extension'. In Section \ref{class2} we will determine the fibre sizes of all word maps over $p$-groups of class $2$ with exponent $p$. We will then give some examples in Section \ref{examples}.

\section{Enumerating sizes of fibres}\label{mainsec}

Let $G$ be a finite $p$-group and, for each $i\in\mathbb{N}$, write $\Irr^i(G) = \{$irreducible complex characters of $G$ of degree $p^i\}$ and $\Irr(G)$ for the set of all irreducible complex characters.
For a complex variable $s$ and $g\in G$ write
\begin{equation}\label{twistedzeta}
\zeta^i_G(s,g) = \sum_{\chi\in\Irr^i(G)}\frac{\chi(g)}{\chi(1)^s}.
\end{equation}
We will also write
\begin{equation}\label{twistedzetamain}
\zeta_G(s,g) = \sum_{\chi\in\Irr(G)}\frac{\chi(g)}{\chi(1)^s}=\sum_{i\in\mathbb{N}}\zeta^i_G(s,g).
\end{equation}
It is clear from equations (\ref{equation1}) and (\ref{twistedzetamain}) that
$$
N^G_t(g) = |G|^{2t-1}\zeta_G(t,g)
$$
and that
\begin{equation}\label{equationprobzeta}
P^G_t(g) = \frac{1}{|G|}\zeta_G(t,g).
\end{equation}
In this section we prove Theorem \ref{main} which allows us to compute expressions like equations (\ref{twistedzetamain}) and (\ref{equationprobzeta}) for finite $p$-groups of nilpotency class less than $p$ with exponent $p$.

Note that when $s=1$ and $g=1$, the identity of $G$, the sum $\zeta_G(1,1)$ is simply the class number $k(G)$ of $G$ and that if $g$ is not in the derived group of $G$ then $\zeta_G(s,g)=0$ since $N_s^G(g)=0$, see equation (\ref{equation1}). An analogue of equation (\ref{twistedzetamain}), a \textit{twisted} zeta function, is studied by Jaikin-Zapirain c.f.\ \cite[Theorem 1.2]{Jaikinzeta} where he shows that it is a rational function in $p^{-s}$ for any $g\in G$ where $G$ is a FAb uniform pro-$p$ group.

We also note that the twisted zeta function in equation (\ref{twistedzetamain}) is multiplicative under direct products in the following sense: for $g_1\in G_1$ and $g_2\in G_2$ where $G_1$ and $G_2$ are groups we have
$$
\zeta_{G_1\times G_2}(s,g_1 g_2) = \zeta_{G_1}(s,g_1).\zeta_{G_2}(s,g_2).
$$

In \cite{VollOBrien} O'Brien \& Voll use the Kirillov orbit method to enumerate the irreducible complex characters, of each degree, of finite $p$-groups of nilpotency class less than $p$. If the group $G$ is of exponent $p$, then the number of characters, of each degree, of $G$ can be described in terms of the number of rational points of certain algebraic varieties.
We can analogously compute the sums $\zeta_G^i(s,g)$ by counting the numbers of rational points of the algebraic varieties considered by O'Brien \& Voll that intersect a hyperplane characterized by $g$.

We now fix a $p$-group $G$ of nilpotency class less than $p$ with exponent $p$ and an element $g\in G'$.

The Lazard correspondence establishes an order-preserving bijection between finite $p$-groups of nilpotency class $c<p$ and finite nilpotent Lie rings of $p$-power order and class $c<p$; cf.\ \cite[Example 10.24]{Khukhro}. Let $c<p$ be the nilpotency class of $G$. Let $\mathfrak{g}=\log(G)$ be the finite Lie ring associated to $G$ by the Lazard correspondence. The Kirillov orbit method gives a correspondence between characters of $G$ and orbits in $\hat{\mathfrak{g}}:=\Hom_{\mathbb{Z}}(\mathfrak{g},\mathbb{C}^*)$, the Pontryagin dual of $\mathfrak{g}$, under the co-adjoint action of $G$ on $\hat{\mathfrak{g}}$; cf.\ \cite[Theorem 2.6]{BoySab} or \cite[Theorem 4.4]{Jon}. Under this correspondence each orbit $\Omega$ of size, say, $p^{2i}$ gives rise to a character $\chi_{\Omega}$ of degree $p^i$ and all characters are of this form. We hav
$$
\zeta^i_G(s,g)=\sum_{\Omega\subseteq\hat{\mathfrak{g}}, |\Omega|=p^{2i}}\frac{\chi_{\Omega}(g)}{\chi_{\Omega}(1)^s},\\
$$
where we sum over orbits $\Omega$ of $\hat{\mathfrak{g}}$ and $\chi_{\Omega}$ is the character of $G$ that corresponds to the orbit $\Omega$. For each $\omega\in\Omega$ we denote $\chi_{\omega} = \chi_{\Omega}$.

For each $\omega\in\hat{\mathfrak{g}}$ we write $B_{\omega}$ for the bi-additive, skew-symmetric form $B_{\omega}:\mathfrak{g}\times\mathfrak{g}\rightarrow\mathbb{C}^*$, $(u,v)\mapsto\omega([u,v])$ and $\Rad(B_{\omega})$ for the radical of $B_{\omega}$. Note that the centre $\mathfrak{z}$ of $\mathfrak{g}$ is contained in $\Rad(B_{\omega})$ and the form $B_{\omega}$ only depends on the restriction of $\omega$ to $\mathfrak{g}'$. From \cite{VollOBrien} we have, for each $i$,
\begin{eqnarray*}
\zeta^i_G(s,g)&=&p^{-2i}\sum_{\omega\in\hat{\mathfrak{g}}, |\mathfrak{g}:\Rad(B_{\omega})|=p^{2i}}\frac{\chi_{\omega}(g)}{\chi_{\omega}(1)^s}\\
&=&|\mathfrak{g}/\mathfrak{g'}|p^{-2i}\sum_{\omega\in\widehat{\mathfrak{g'}}, |\mathfrak{g}:\Rad(B_{\omega})|=p^{2i}}\frac{\chi_{\omega}(g)}{\chi_{\omega}(1)^s}\\
&=&|\mathfrak{g}/\mathfrak{g'}|p^{-2i}\sum_{\omega\in\widehat{\mathfrak{g'}}, |\Rad(B_{\omega}):\mathfrak{z}|=p^{-2i}|\mathfrak{g}/\mathfrak{z}|}\frac{\chi_{\omega}(g)}{\chi_{\omega}(1)^s}.\\
\end{eqnarray*}
By the Kirillov orbit method we have, for all $\omega\in\hat{\mathfrak{g}}$ such that $|\mathfrak{g}:\Rad(B_{\omega})|=p^{2i}$,
$$
\chi_{\omega}(g) = p^{-i}\sum_{\upsilon\in\Omega_{\omega}}\upsilon(g),
$$
where $\Omega_{\omega}$ is the orbit containing $\omega$ and we identify $g$ with an element of $\mathfrak{g}=\log(G)$. Hence
$$
\zeta^i_G(s,g)=|\mathfrak{g}/\mathfrak{g'}|p^{-i(3+s)}\sum_{\omega\in\widehat{\mathfrak{g'}}, |\Rad(B_{\omega}):\mathfrak{z}|=p^{-2i}|\mathfrak{g}/\mathfrak{z}|}\sum_{\upsilon\in\Omega_{\omega}}\upsilon(g).
$$

Theorem B in \cite{VollOBrien} gives a geometric characterization of this sum in terms of the numbers of certain rational points of rank varieties of matrices of linear forms. Assume that $\mathfrak{o}$ is a compact discrete valuation ring of characteristic zero with maximal ideal $\mathfrak{p}$ and residue field $\bold{k} = \mathfrak{o}/\mathfrak{p}$ of characteristic $p$. Suppose that $\mathfrak{g}$ is a finite, nilpotent $\mathfrak{o}$-Lie algebra of class $c<p$ and that both $\mathfrak{g}/\mathfrak{z}$ and $\mathfrak{g}'$ are annihilated by $\mathfrak{p}$. Set $a:=\rk_{\bold{k}}(\mathfrak{g}/\mathfrak{z})$ and $b:=\rk_{\bold{k}}(\mathfrak{g'})$ and fix bases $\bold{e}=\{e_1,...,e_a\}$ and $\bold{f}=\{f_1,...,f_b\}$ for $\mathfrak{g}/\mathfrak{z}$ and $\mathfrak{g'}$ respectively. Choose structure constants $\lambda_{ij}^k\in\bold{k}$ such that
$$
[e_i,e_j] = \sum_{k=1}^{b}\lambda_{ij}^kf_k
$$
and $\lambda_{ij}^k = -\lambda_{ji}^k$ for all $i,j\in\{1,...,a\}$ and $k\in\{1,...,b\}$. Let $\bold{Y} = (Y_1,...,Y_b)$ be independent variables and define the `commutator matrix' (with respect to $\bold{e}$ and $\bold{f}$) $B(\bold{Y})\in\Mat_a(\bold{k}[\bold{Y}])$ given by
$$
B(\bold{Y})_{ij}:=\sum_{k=1}^{b}\lambda_{ij}^{k}Y_k
$$
for all $i,j\in\{1,...,a\}$. For $\bold{y}=(y_1,...,y_b)\in\bold{K}^b$, where $\bold{K}$ is any finite extension of $\bold{k}$, write $B(\bold{y})\in\Mat_a(\bold{K})$ for the matrix obtained by evaluating the variables $Y_i$ at $y_i$. Note that the matrices $B(\bold{y})$ are skew-symmetric, have even rank
and that $\det(B(\bold{Y}))$ is a square in $\bold{k}[\bold{Y}]$, whose square root $\Pf(B(\bold{Y})):=\sqrt{\det(B(\bold{Y}))}$ is the \textit{Pfaffian} of $B(\bold{Y})$. If $a$ is odd then $\Pf(B(\bold{Y}))=0$.

Fix a non-trivial additive character $\phi:\bold{k}\rightarrow\mathbb{C}^*$. For $a\in\bold{k}$ define $\phi_a(x) = \phi(ax)$. The map $a\mapsto\phi_a$ is an isomorphism between $\bold{k}$ and its Pontryagin dual $\widehat{\bold{k}}$.
We have an isomorphism between $\mathfrak{g}'$ and its dual $\widehat{\mathfrak{g}'}$ and also a canonical isomorphism between $\mathfrak{g}'$ and its linear dual $\Hom_{\mathbb{F}_q}(\mathfrak{g}',\bold{k})$. We now fix an isomorphsim $\psi_1:\widehat{\mathfrak{g}'}\rightarrow\Hom_{\bold{k}}(\mathfrak{g}',\bold{k})$. The dual $\bold{k}$-basis $f^{\vee} = (f_k^{\vee})$ for $\Hom_{\bold{k}}(\mathfrak{g}',\bold{k})$ gives a coordinate system
\begin{eqnarray*}
\psi_2:\Hom_{\bold{k}}(\mathfrak{g}',\bold{k})&\rightarrow&\bold{k}^b;\\
y=\sum_{k=1}^{b}y_kf_f^{\vee}&\mapsto&\bold{y}=(y_1,...,y_b).
\end{eqnarray*}
Set $\psi:=\psi_2\circ\psi_1:\widehat{\mathfrak{g}'}\rightarrow\bold{k}^b$, an isomorphism. Under this isomorphism we may identify elements $\bold{y}\in\bold{k}^b$ with elements $\omega_{\bold{y}}\in\widehat{\mathfrak{g}'}$.

Assume that $\mathfrak{O}$ is an unramified extension of $\mathfrak{o}$, with maximal ideal $\mathfrak{P}$. Write $\mathfrak{g}(\mathfrak{O})$ for $\mathfrak{g}\otimes_{\mathfrak{o}}\mathfrak{O}$ and $\mathfrak{z}(\mathfrak{O})$ for $\mathfrak{z}\otimes_{\mathfrak{o}}\mathfrak{O}$. We identify the residue field $\mathfrak{O}/\mathfrak{P}$, a finite extension of $\bold{k}$, with $\mathbb{F}_q$. The derived $\mathfrak{O}$-Lie algebra $\mathfrak{g}(\mathfrak{O})'$ as well as the quotient $\mathfrak{g}(\mathfrak{O})/\mathfrak{z}(\mathfrak{O})$ are annihilated by $\mathfrak{P}$. We denote the corresponding bases, $\bold{e}\otimes_{\bold{k}}1$ for $\mathfrak{g}(\mathfrak{O})/\mathfrak{z}(\mathfrak{O})$ and $\bold{f}\otimes_{\bold{k}}1$ for $\mathfrak{g}(\mathfrak{O})'$, obtained by this base extension by $\bold{e}$ and $\bold{f}$ respectively. We consider $\bold{e}$ and $\bold{f}$ as $\mathbb{F}_q$ bases for the respective $\mathbb{F}_q$-vector spaces of dimensions $a$ and $b$. The group $G=G(\mathfrak{O}):=\exp(\mathfrak{g}(\mathfrak{O}))$ has the property that both $G(\mathfrak{O})'$ and $G(\mathfrak{O})/Z(G(\mathfrak{O}))$ have exponent $p$. It follows from \cite[Theorem B]{VollOBrien} that
$$
\zeta^i_G(s,g)=|\mathfrak{g}/\mathfrak{g'}|p^{-i(3+s)}\sum_{\bold{y}\in\mathbb{F}_q^b,  \rk(B(\bold{y}))=2i}\sum_{\bold{v}\in\Omega_{\omega_\bold{y}}}\bold{v}(g).
$$
Since $|\Omega_{\omega_\bold{y}}| = p^{2i}$ we have
$$
\zeta^i_G(s,g)=|\mathfrak{g}/\mathfrak{g'}|p^{-i(1+s)}\sum_{\bold{y}\in\mathbb{F}_q^b,  \rk(B(\bold{y}))=2i}\omega_{\bold{y}}(g).
$$
Let
\begin{equation}\label{eqnK}
K^i_G(g)=\#\{\bold{y}\in\mathbb{F}_q^b:  \rk(B(\bold{y}))=2i, g\in\Ker(\omega_{\bold{y}})\}
\end{equation}
and
\begin{equation}\label{eqnV}
V^i_G(g)=\#\{\bold{y}\in\mathbb{F}_q^b:  \rk(B(\bold{y}))=2i, \ord(\omega_{\bold{y}}(g))\neq 1\}.
\end{equation}
Note that $\sum_i K^i_G(1) = q^b$ and $V_G^i(1)=0$ for all $i$ whilst for $1\neq g\in G$ we have $\sum_i K^i_G(g)=q^{b-1}$.

We write $K_G(g)$ and $V_G(g)$ for the vectors $(K^i_G(g))_i$ and $(V^i_G(g))_i$ respectively. The numbers $\omega_{\bold{y}}(g)$ are $p$-th roots of unity. Since the sum of the $q-1$ Galois conjugates of non-trivial $q$-th roots of unity is $-1$ we have the following:
\begin{theorem}\label{main}
Let $\mathfrak{o}$ be a compact discrete valuation ring of characteristic zero with residue field $\bold{k}$ of characteristic $p$ and let $\mathfrak{g}$ be a finite, nilpotent $\mathfrak{o}$-Lie algebra of class $c<p$. Assume that $\mathfrak{g}'\cong\bold{k}^b$ and $\mathfrak{g}/\mathfrak{z}\cong\bold{k}^a$ as $\bold{k}$-vector spaces. Let $\mathfrak{O}$ be a finite, unramified extension of $\mathfrak{o}$, with residue field isomorphic to $\mathbb{F}_q$. Write $G:=G(\mathfrak{O})$ for the group associated with $\mathfrak{g}(\mathfrak{O})$ under the Lazard correspondence. Then, for each $i$,
$$
\zeta^i_G(s,g)=|G/G'|p^{-i(1+s)}(K^i_G(g)-\frac{1}{q-1}V^i_G(g)),
$$
where $K^i_G$ and $V^i_G$ are defined in equations (\ref{eqnK}) and (\ref{eqnV}).
\end{theorem}
\noindent This follows immediately from the discussion above.

\section{Nilpotent class $2$ groups}\label{class2}

Let $G$ be a finite group, $w(x_1,...,x_n)$ a group word and $w$ the associated word map. Recall that $N_w^G(g)$ is the number of solutions to $w= g$
and that $P_w^G(g)$ is the probability that a random $n$-tuple $\textbf{g}=(g_{1},...,g_{n})\in G^{(n)}$ satisfies $w(\textbf{g})=g$.
The following corollary follows immediately from the results in \cite{Levy1}:

\begin{corollary}[\cite{Levy1}]
Let $G$ be a finite $p$-group of nilpotency class $2$ with exponent $p$ and let $w$ be a group word. Then there exists a word $v$ of the following form:
\begin{itemize}
\item[i)] $v = x$; or
\item[ii)] $v = c_t$ for some $t$,
\end{itemize}
such that $P_w^G(g) = P_v^G(g)$ for all $g\in G$.
\end{corollary}

By Theorem \ref{main} we have determined the fibre sizes of all word maps over all $p$-groups of nilpotency class $2$ with exponent $p$:

\begin{corollary}
Let $p$ be an odd prime and let $G$ be a finite $p$-group of nilpotency class $2$ with exponent $p$ and let $w$ be a group word. Then either
\begin{itemize}
\item[i)] $P_w^G(g) = \frac{1}{|G|}$ for all $g\in G$; or
\item[ii)] $P_w^G(g) = \frac{1}{|G|}\zeta_G(t,g)$ for some $t$ depending on $w$.
\end{itemize}
\end{corollary}

\section{Examples}\label{examples}

We compute the sums $\zeta^i_G(s,g)$ for various relatively free $p$-groups with exponent $p$. Note that for $i=0$, the only vector $\bold{y}$ giving rise to a matrix $B(\bold{y})$ such that $\rk(B(\bold{y}))=0$ (see equations (\ref{eqnK}) and (\ref{eqnV})) is $\bold{y} = \bold{0}$. Suppose that $i\neq 0$. Since the rank of the matrix $B(\bold{y})$, for given $\bold{y}\in\mathbb{F}_q^b$, is invariant under scalar multiplication the rank of the matrix $B(\tilde{\bold{y}})$ with $\tilde{\bold{y}}=(\tilde{y_1}:...:\tilde{y_b})\in\mathbb{P}^{b-1}(\mathbb{F}_q)$ is well defined. Under the Lazard correspondence we may identify $1\neq g\in G'$ with an element $\bold{g}\in\mathfrak{g}'\cong\mathbb{F}_q^{b}$, the associated Lie algebra, and similarly identify $\bold{g}$ with $\tilde{\bold{g}}\in\mathbb{P}^{b-1}(\mathbb{F}_q)$. Since the matrix $B$ is skew-symmetric its determinant $\det(B)$ is a square whose square root $\Pf(B):=\sqrt{\det(B)}$ is the \textit{Pfaffian} of $B$. If $a$ is odd, then $\Pf(B) = 0$. Assume that $\Pf(B) \neq 0$. Then $\Pf(B)$ defines a hypersurface in $\mathbb{P}^{b-1}$ and the $\mathbb{F}_q$-rational points $\tilde{\bold{y}}$ of this hypersurface correspond to matrices $B(\tilde{\bold{y}})$ of a certain non-maximal rank. The $\mathbb{F}_q$-rational points which do not lie on the hypersurface will be of maximal rank $       a$. We refer to points $\tilde{\bold{y}}\in\mathbb{P}^{b-1}(\mathbb{F}_q)$ as being `of rank 2i' for some $i$ if the associated matrix $B(\tilde{\bold{y}})$ is of rank $2i$. The condition $g\in\Ker(\omega_{\bold{y}})$ in the definition of $K^i_G(g)$ (see equation (\ref{eqnK})) defines a hyperplane $H_g$ given by the dot product $\tilde{\bold{g}}.\tilde{\bold{y}}=0$ in $\mathbb{P}^{b-1}(\mathbb{F}_q)$. We may thus talk about the `hyperplane defined by $g$' as $H_g$. The numbers $K^i_G(g)$ are simply the number of $\mathbb{F}_q$-rational points of this hyperplane which intersect the hypersurface defined by $\Pf(B)$ in $\mathbb{P}^{b-1}(\mathbb{F}_q)$ giving rise to a point `of rank $2i$'.

We will now proceed with some examples. In each case it turns out that the matrices $B(\bold{Y})$ have the form
$$
B(\bold{Y}) =
\begin{pmatrix}
0&U(\bold{Y})\\
-U(\bold{Y})^{\tr}& 0
\end{pmatrix}
$$
where $U(\bold{Y})$ is a matrix.

\begin{example}[Heisenberg group, H($\mathbb{F}_q$)]
Suppose that $G$ is the Heisenberg group H($\mathbb{F}_q$). The matrix $U(\bold{Y})$, where $\bold{Y}=(Y_1)$ has a single variable, is simply the $1\times 1$ matrix with entry $Y_1$. There are two cases, $i=0$ and $i=1$, for the rank of the matrix $U(\bold{Y})$. The number of vectors $\bold{y}$ in each case is $1$ and $q-1$ respectively. When $g=1$ is the identity, we have $K^0_G(1)=1$ and $V^0_G(1) = 0$ so that $\zeta^0_G(s,1)=q^2$. We also have $K^1_G(1)=q-1$ and $V^1_G(1) = 0$ so that $\zeta^1_G(s,1)=q^{1-s}(q-1)$. Together, this gives us
$$
\zeta_G(s,1)=q^2+q^{-s+1}(q-1)
$$
and when $s=1$ this is simply $k(G)$, the class number.

Suppose now that $1\neq g\in G'$. This occurs with multiplicity $(q-1)$. It is not hard to see that $\zeta^0_G(s,g)=q^2$ since $K^0_G(g)=1$ and $V^0_G(g) = 0$. Also, $K^1_G(g)=0$ and $V^1_G(g) = q-1$ so that
$$
\zeta_G(s,g)=q^2-q^{-s+1}.
$$

\end{example}

\begin{example}[A quadric surface in $\mathbb{P}^2(\mathbb{F}_q)$]
Let $\mathfrak{g}$ be the $7$-dimensional nilpotent $\mathbb{F}_q$-Lie algebra of class $2$ with $\mathbb{F}_q$-basis $(x_1,...,x_4,y_1,y_2,y_3)$ subject to the relations $[x_1,x_3]=y_1$, $[x_1,x_4]=y_2$, $[x_2,x_3]=y_3$, $[x_2,x_4]=y_1$. With respect to this basis we have
$$
U(\bold{Y}) =
\begin{pmatrix}
Y_1&Y_2\\
Y_3& Y_1
\end{pmatrix}.
$$
The determinant of this matrix is defines the quadric surface $Y_1^2-Y_2Y_3$ in $\mathbb{P}^2(\mathbb{F}_q)$. There are three cases, $i=0,1$ and $2$, for the rank of the matrix $U(\bold{Y})$. The number of such vectors $\bold{y}$ in each case is $1$, $(q+1)(q-1)$ and $q^2(q-1)$ respectively. As usual, the $i=0$ case corresponds to when $\bold{y}$ is zero, $i=1$ when $\tilde{\bold{y}}$ lies on the curve and $i=2$ otherwise. When $g$ is the identity we can compute
$$
K_G(1)=(1,(q+1)(q-1),q^2(q-1))
$$
and
$V_G(1)=(0,0,0)$ so that $k(G)=q^4+q^2(q+1)(q-1)+q^2(q-1)$.

Now suppose that $g\neq 1$. We have two cases. The first case is when $H_g$ corresponds to a tangent of the surface $Y_1^2=Y_2Y_3$. This occurs with multiplicity $(q+1)(q-1)$. In this case
$$
K_G(g)=(1,(q-1),q(q-1))
$$
and
$$
V_G(g)=(0,q(q-1),(q^2-q)(q-1)).
$$

The second case is when $H_g$ corresponds to a line which intersects the surface $Y_1^2=Y_2Y_3$ at two distinct points. This occurs with multiplicity $q^2(q-1)$. In this case
$$
K_G(g)=(1,2(q-1),(q-1)(q-1))
$$
and
$$
V_G(g)=(0,(q-1)(q-1),(q^2-q+1)(q-1)).
$$
\end{example}

\begin{example}[A quadric surface in $\mathbb{P}^3(\mathbb{F}_q)$]
Let $\mathfrak{g}$ be the $8$-dimensional nilpotent $\mathbb{F}_q$-Lie algebra of class $2$ with $\mathbb{F}_q$-basis $(x_1,...,x_4,y_1,...,y_4)$ subject to the relations $[x_1,x_3]=y_1$, $[x_1,x_4]=y_2$, $[x_2,x_3]=y_3$, $[x_2,x_4]=y_4$. With respect to this basis we have
$$
U(\bold{Y}) =
\begin{pmatrix}
Y_1&Y_2\\
Y_3& Y_4
\end{pmatrix}.
$$
The determinant of this matrix defines the quadric surface $Y_1Y_4-Y_2Y_3$ in $\mathbb{P}^3(\mathbb{F}_q)$. There are three cases, $i=0,1$ and $2$, for the rank of the matrix $U(\bold{Y})$. The number of vectors $\bold{y}$ in each case is $1$, $(q+1)^2(q-1)$ and $(q^3-q)(q-1)$ respectively. As usual, the $i=0$ case corresponds to when $\bold{y}$ is zero, $i=1$ when $\tilde{\bold{y}}$ lies on the curve and $i=2$ otherwise. When $g$ is the identity we compute
$$
K_G(1)=(1,(q+1)^2(q-1),q^4-1-(q+1)^2(q-1))
$$
and
$V_G(1)=(0,0,0)$ so that $k(G)=q^4+q^2(q+1)^2(q-1)+q^4-1-(q+1)^2(q-1)$.

Now suppose that $g\neq 1$. We have two cases. The first case is when $g$ corresponds to a tangent of the surface $Y_1Y_4=Y_2Y_3$. This occurs with multiplicity $(q+1)^2(q-1)$. In this case
$$
K_G(g)=(1,(2q+1)(q-1),(q^2-q)(q-1))
$$
and
$$
V_G(g)=(0,q^2(q-1),(q^3-q^2)(q-1)).
$$

The second case is when $g$ corresponds to a plane which intersects the surface $Y_1Y_4=Y_2Y_3$ and is non-tangent. This occurs with multiplicity $(q^3-q)(q-1)$. In this case
$$
K_G(g)=(1,(q+1)(q-1),q^2(q-1))
$$
and
$$
V_G(g)=(0,q(q+1)(q-1),(q^3-q^2-q)(q-1)).
$$
\end{example}

\noindent The following examples were studied by Boston \& Isaacs \cite{BostonIsaacs}.

\begin{example}[Elliptic curves in $\mathbb{P}^2(\mathbb{F}_p)$]
Let $p$ be a prime and $\alpha\in\mathbb{F}_p^*$. Let $\mathfrak{g}_{\alpha}$ be the $9$-dimensional nilpotent $\mathbb{F}_p$-Lie algebra of class $2$ with $\mathbb{F}_p$-basis $(x_1,...,x_6,y_1,y_2,y_3)$ subject to the relations $[x_1,x_4]=y_1$, $[x_1,x_5]=y_2$, $[x_1,x_6]=\alpha y_3$, $[x_2,x_4]=y_3$, $[x_2,x_5]=y_1$, $[x_2,x_6]=y_2$, $[x_3,x_4]=y_3$, $[x_3,x_6]=y_1$. With respect to this basis we have
$$
U(\bold{Y}) =
\begin{pmatrix}
Y_1&Y_2&\alpha Y_3\\
Y_3& Y_1&Y_2\\
Y_3&0&Y_1
\end{pmatrix}.
$$
The determinant of this matrix defines an elliptic curve $E_{\alpha}$ in $\mathbb{P}^2(\mathbb{F}_p)$ and let $n_{\alpha}$ denote the number of $\mathbb{F}_p$-rational points of the curve $E_{\alpha}$. There are three cases, $i=0,1$ and $2$, for the rank of the matrix $U(\bold{Y})$. The number of such vectors $\bold{y}$ in each case is $1$, $n_{\alpha}(p-1)$ and $(p^2+p+1-n_{\alpha})(p-1)$ respectively. As usual, the $i=0$ case corresponds to when $\bold{y}$ is zero, $i=1$ when $\tilde{\bold{y}}$ lies on the curve and $i=2$ otherwise. When $g$ is the identity we can compute
$$
K_G(1)=(1,n_{\alpha}(p-1),(p^2+p+1-n_{\alpha})(p-1))
$$
and $V_G(1)=(0,0,0)$ so that $k(G)=p^6+p^2n_{\alpha}(p-1)+(p^2+p+1-n_{\alpha})(p-1)$.

Now suppose that $g\neq 1$. Let $k_{\alpha}$ denote the number of inflection points of $E_{\alpha}$. We have four cases depending on how many times the line $H_g$ corresponding to $g$ intersects the elliptic curve $E_{\alpha}$ at rational points, this can be zero, once, twice or three times occurring with multiplicities $p^2+p+1-n_{\alpha}(p+1)+\frac{2}{3}{n_{\alpha} \choose 2}+\frac{n_{\alpha}-k}{3}$, $n_{\alpha}(p+1)-(n_{\alpha}-k)-{n_{\alpha} \choose 2}$, $(n_{\alpha}-k)(p-1)$ and ${n_{\alpha} \choose 2}$ respectively.

Suppose that $1\neq g$ defines a line which intersects the curve at $m$ rational points where $m=0,1,2,3$. We have
$$
K_G(g)=(1,m(p-1), (p+1-m)(p-1))
$$
and
$$
V_G(g)=(0,(n_{\alpha}-m)(p-1),(p^2+m-n_{\alpha})(p-1)).
$$
\end{example}

\section{Acknowledgements}
I would like to thank Christopher Voll for his support, guidance and insight throughout this project and for numerous helpful discussions.

\bibliographystyle{plain}


\section{introduction}
Graphene\cite{1,2} exhibits remarkable features such as massless
Dirac fermions, an anomalous Berry's phase, and signatures of
Hofstadter's spectrum,\cite{3,4,5,6} while allowing for the gate-induced carrier type conversion, from
holes to electrons, without crossing a
bandgap.\cite{7,8,9,10} When the Fermi level is placed at
the Dirac point in monolayer graphene, one expects the
\textit{density} of carriers, $n_{q}$, i.e., electrons ($n_e$) and
holes ($n_h$), to vanish, leading to a divergent Hall resistance,
$R_{xy}$, and a diverging resistivity $\rho_{xx}$, since $R_{xy}$
and $\rho_{xx}$ are  $\propto 1/n_{q}$. Remarkably, experiment has
reported a finite, nearly quantized $\rho = \sigma^{-1} \approx
h/4e^{2}$, reminiscent of the prediction $\sigma_{min} = e^2/h$
per massless Dirac channel,\cite{7,9} electron-hole
puddles,\cite{9,12b,14,13} along with Hall effect compensation.\cite{1} Here, we examine vanishing- instead of
diverging- Hall effect over the p$\leftrightarrow$n transition, in CVD
graphene,\cite{15} in addition to $R_{xy} \rightarrow 0$, as  the diagonal resistance through the Dirac peak indicates
$\rho_{xx} \approx h/4e^{2}$ at the nominal neutrality point. We
reproduce these characteristics using a parabolic distribution $f(V_{N})$ of
neutrality potentials, $V_{N}$, in an ambipolar multi-conduction
model.\cite{16}  The results serve to extract the effective
residual electron- and hole- densities at the nominal neutrality point, and as a function of the gate voltage around the neutrality point. The observed simultaneous non-vanishing electron- and hole- densities in this zero bandgap material at liquid helium temperatures suggest that the applied current could be simultaneously
carried by streams of electrons and holes within spatially separated current domains in these large area CVD graphene devices.

\begin{figure}[t]
\begin{center}
\includegraphics[width=60mm]{ordhallgraphene_fig1v2eps}
\end{center}
\caption{(Color online) (a) The Raman intensity vs. the Raman shift for
graphene-on-$SiO_{2}/Si$-substrate. The D-, G-, and 2D- peaks have
been marked. The 2D-to-G intensity ratio, $I_{2D}/I_{G} = 2$,
suggests monolayer graphene. (b) The diagonal resistance,
$R_{xx}$, and the Hall resistance, $R_{xy}$, have been plotted vs.
the magnetic field, $B$, for a gate-voltage $V_{G} = 73$ V. The
figure indicates a relatively large ($ \approx 8$ per-cent) weak
localization correction. \label{fig: epsart}}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=50mm]{ordhallgraphene_fig2v2eps.eps}
\end{center}
\caption{ (Color online) (a) At $T=1.5$K, $R_{xx}$ is plotted vs. $B$ with $V_{G}$ as parameter. Here,
$R_{xx}$ initially
increases up to $V_{G} = 75 $ V. For $V_{G} \ge 80$ V, $R_{xx}$
decreases with increasing $V_{G}$. (b) Fits to graphene weak
localization theory, see text.
(c) The fit extracted zero-field resistance, $R_{xx}^{0}$, and the
amplitude of the weak localization correction, $\Delta R_{xx}$,
are shown vs. $V_{G}$. (d) The fit extracted phase coherence
field, $B_{\phi}$ and the phase coherence length $l_{\phi}$ are
exhibited vs. $V_{G}$.  (e) The fit extracted inter-valley field,
$B_{iv}$ and the inter-valley length $l_{iv}$ are exhibited vs.
$V_{G}$. (f) At $T=1.5$K, the Hall resistance $R_{xy}$ is plotted
vs. $B$ with $V_{G}$ as a parameter. (g) Note the non-monotonic
variation in $dR_{xy}/dB$ vs. $V_{G}$.
\label{afig: epsart}}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[width=60mm]{ordhallgraphene_fig3eps}
\end{center}
\caption{ (Color online) (a) The Hall resistance, $R_{xy}$, is exhibited vs. $B$
for different $V_{G}$. The green (red) traces show a positive slope
$dR_{xy}/dB$ that increases (decreases) with $V_{G}$. The cyan
trace marks the boundary between the red and dark
blue traces. Dark blue traces show a negative slope. (b) $dR_{xy}/dB$ evaluated from
panel (a) are plotted as red symbols vs. $V_{G}$.  The green and
brown symbols show $1/((dR_{xy}/dB)e)$, the
carrier density in a single carrier model. (c) $R_{xx}$ (blue trace) peaks at the nominal neutrality
point (dashed vertical line). The evaluated
$\sigma_{xy}/(\sigma_{xx} B)$ (green symbols), would suggest a drop in the hole mobility over the range
$55 \le V_{G} \le 72$V.
\label{figb:epsart}}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[width=60mm]{ordhallgraphenefig4_fig4v4eps}
\end{center}
\caption{ (Color online) (a) The blue trace shows the distribution function,
$f(V_{N})$, of the neutrality voltages. The red and magenta traces show the effective
hole and electron densities, respectively. The black
lines show the hole and electron densities in ideal
monolayer graphene. (b) The data (blue symbols) and best fit (brown
curve) of the $dR_{xy}/dB$ vs. $V_{G}$ obtained for the effective
carrier densities of panel (a). Inset: A comparison of
the $V_{G}$ dependence of $dR_{xy}/dB$ in ideal monolayer
graphene and CVD graphene. (c) The data (green symbols) and fit
(red curve) of the $R_{xx}$ vs. $V_{G}$ obtained for the effective
carrier densities of panel (a). Inset: Near the nominal
neutrality point, coexisting electron and hole streams, with
oppositely oriented Hall electric fields ($E_{H,e}$ and
$E_{H,h}$), carry the applied current. The Hall effect then
vanishes due to voltage compensation. \label{bfig: epsart}}
\end{figure}
Monolayer graphene was prepared using the CVD technique,\cite{15}
processed into mm-scale devices, and measured in a low temperature
magnet cryostat, see supplementary material.\cite{15b} The inset of Fig. 1(a) shows the graphene on $SiO_{2}/Si$, while
figure 1(a) exhibits the Raman scattering spectrum. Here, the
D-peak is small and the 2D-to-G intensity ratio, $I_{2D}/I_{G} =
2$, suggests monolayer graphene.\cite{15} Figure 1(b) shows the
diagonal resistance, $R_{xx}$, and $R_{xy}$ vs. the magnetic
field, $B$, at a back gate-voltage $V_{G} = 73$ V, close to the
nominal neutrality voltage, $\nu = 72.8 V$.
This panel (Fig. 1(b)) shows a large ($ \approx 8\%$) weak
localization correction to
$R_{xx}$.\cite{20,8,21,21b,22,9,10} The $|R_{xy}| = 3.7$
ohms at $B=0.5$ Tesla evaluates to $B/R_{xy}e = 8.5 \times
10^{13}$ cm$^{-2}$ at a $V_{G}$ where one expects $n_{e} < 1\times
10^{12}cm^{-2}$.

The next figures provide an experimental overview of the approach
to a vanishing Hall effect state at $V_G=\nu$. Figure 2(a) shows that $R_{xx}$ changes non-monotonically with
$V_{G}$, as $R_{xx}$ initially increases, followed by a decrease
above $V_{G} = 80 $ V. Here, $\nu \approx 77$ V. Also observable
in Fig. 2(a) is the weak localization correction which increases
towards $V_G = \nu$. Fig. 2(b) exhibits good data fits to weak
localization theory,\cite{21b} see supplementary \cite{15b}, on both the nominally hole- and nominally
electron-side  of $V_{G} = \nu$.  Fig.
2(c) shows the fit extracted zero-field resistance, $R_{xx}^{0}$,
and the weak localization amplitude, $\Delta R_{xx}$, vs. $V_{G}$.
At $V_{G}=75V$, $R_{xx}^{0} = 7.1k\Omega$ and $R_{xx} (0.4 T) =
6.5 k\Omega$ lie within $\approx 10 \%$ of $h/4e^{2} = 6.45
k\Omega$. Fig. 2(d) exhibits the fit extracted phase coherence
field $B_{\phi}$ and the phase coherence length $l_{\phi}$ vs.
$V_{G}$. The fit extracted inter-valley
field, $B_{iv}$ and the inter-valley length $l_{iv}$ are exhibited
vs. $V_{G}$ in Fig. 2(e). In sum, the fits work equally well on
either side of the nominal Dirac point, and the fit parameters,
exhibited in Fig. 2(c)-2(e), suggest increasing "localization"
with gate-induced carrier depletion.

Fig. 2(f) reports the Hall resistance. Here, the slope of $R_{xy}$ vs. $B$ is positive and
increases with $V_{G}$ up to $V_{G} = 60V$. From $60 \le V_{G} \le
75 V$, the slope remains positive although now it decreases in
magnitude with increasing $V_{G}$. Above $V_{G} = 75V$, the slope
turns negative and increases in magnitude up to $V_{G} = 100 V$.
Finally, the negative slope decreases in magnitude above $V_{G}
=100 V$. Fig. 2(g) shows the non-monotonic variation in
$dR_{xy}/dB$ vs. $V_{G}$, where $dR_{xy}/dB = \Delta R_{xy}/\Delta
B = R_{xy} (B)/B$ since $R_{xy}(0)
= 0$. Here, the dashed lines mark an anomalous band where the
$|dR_{xy}/dB|$ decreases in magnitude towards $V_G= \nu$ although
$n_q$ should decrease in this direction.

Characteristics over the crossover of the nominal neutrality point are
examined at smaller $V_{G}$ intervals in Fig. 3. In Fig. 3(a),
for the green traces, the positive slope $dR_{xy}/dB$ increases
with $V_{G}$. For the red traces, the positive slope decreases
with increasing $V_{G}$. The cyan trace, $V_{G} = 73 V$ marks a
boundary between the red and dark blue traces. [Data at this
$V_{G}$ appear in Fig. 1(b)]. For the dark blue traces, the slope
is negative and it increases in magnitude with increasing $V_{G}$.
$dR_{xy}/dB$ evaluated from Fig. 3(a) are plotted as red symbols
vs. $V_{G}$ in Fig. 3(b). Magenta symbols represent other data not
shown. Note that the $dR_{xy}/dB$ increases with increasing
$V_{G}$ only up to $V_{G} \approx 60V$. Above $V_{G} = 60 V$, the
slope decreases with increasing $V_{G}$ as it vanishes slightly
below $V_{G}= 73$ V. For a further increase in $V_{G}$,
$dR_{xy}/dB$ becomes negative and increases in magnitude with
$V_{G}$. Also shown  with the green and brown symbols in Fig. 3(b), is
$1/((dR_{xy}/dB)e)$, where $1/((dR_{xy}/dB)e)$ would equal the
carrier density for single carrier conduction. Here, $1/((dR_{xy}/dB)e)$ decreases up to $V_{G}
\approx 50$V as expected for hole depletion with increasing
positive $V_G$ via $n_{h} \approx C(\nu - V_{G})$ where $C \approx
7.2 \times 10^{-10} cm^{-2}/V$. The trace then deviates from
expectations as $1/((dR_{xy}/dB)e)$ increases with increasing
$V_{G}$ between $60 \le V_{G} \le 72$V, although one expects
further hole depletion. Similarly, from slightly below $V_{G} =
73$V, increase in $V_{G}$ should increase $n_e$, while the
$1/((dR_{xy}/dB)e)$ (brown symbols) indicates opposite behavior.

The $R_{xx}$ vs. $V_{G}$ data at $B=0$ are shown as the blue trace
in Fig. 3(c). The peak $R_{xx} \approx 6.5 k\Omega \approx
h/4e^{2}$ at $V_G = \nu$. We evaluate $ \sigma_{xy}/(\sigma_{xx}
B)$, where $\sigma_{xx}$ and $\sigma_{xy}$ are the conductivity
tensor components, and plot the result vs. $V_G$, see the green trace,
in Fig. 3(c). This plot suggests that holes in this specimen
exhibit an approximately constant mobility over the range $0 \le
V_{G} \le 50V$. Then, above $V_{G} = 50 V$,
$\sigma_{xy}/(\sigma_{xx} B)$ decreases with increasing $V_{G}$ as
it vanishes at $\nu$. At higher $V_{G}$, $\sigma_{xy}/(
\sigma_{xx} B)$ changes sign indicating electron conduction as it
increases in magnitude with $V_{G}$. Thus, a single carrier
interpretation which associates $\sigma_{xy}/( \sigma_{xx} B)$
with $\mu$ would suggest, unphysically, that the mobility changes with $V_G$.

Since a single carrier interpretation leads to un-physically large
$n_q$ close to the Dirac neutrality point, where $n_{q}$ ought
to vanish, and setting $\sigma_{xy}/(\sigma_{xx} B) = \mu$ also
leads to $\mu$ varying un-physically with $V_{G}$, see Fig. 3(c),
we considered, in the first iteration, ambipolar
conduction in a uniform-carrier-density, two carrier Drude model,
where at each $V_G$, there occurs both an electron and a hole
contribution at a constant density to the conductivities $\sigma_{xx}$ and $\sigma_{xy}$, see supplementary.\cite{15b}
Notably, this approach did not help to produce good fits of the data.

The large neutrality voltage, $V_N$, observed in CVD graphene on
$SiO_{2}/Si$ results from the combined gating effects of trapped
charge in the $SiO_2/Si$ substrate, charge donating impurities
below the graphene,\cite{31} and adsorbates on the graphene layer.
It appeared plausible that such disorder effects could combine to produce a
variation in the neutrality voltage across the mm-scale
CVD graphene specimen. Thus, we introduced a distribution
function, $f(V_{N})$, for $V_{N}$ with $\int f(V_{N})dV_{N} =1$, see supplementary.\cite{15b} Here, we focus upon the
results obtained with the parabolic distribution, $f(V_{N}) = A(w^{2} - (V_N - \nu)^2)$, blue curve
of Fig. 4(a), which is non-vanishing only between
$V_{N} = \nu \pm w$.

A distribution function $f(V_{N})$
introduces a gradation in the neutrality voltage, and a non-uniformity in the carrier density across the specimen. Suppose, for the sake of discussion, that the spread in neutrality voltages
occurs all in one direction, say, the x-direction, across the
width (W) of a Hall bar device of length L, with the left edge
of the Hall bar at $x=0$ is at $V_N = \nu -w$ and the right edge
of the Hall bar at $x=W$ is at $V_N = \nu + w$. Suppose the step
size in $V_{N}$ is $\Delta V_{N}$, then the spread in $V_N$ would
correspond to a number of strips of constant $V_N$, which
percolate from the bottom end of the Hall bar at $y=0$ to the top
end at $y=L$. (Scanning probe studies, see for example Fig. 3(a) of ref. \cite{13},
have suggested percolation of the hole and electron puddles over micron length scales). The spatial width, $\Delta x_{i}$, of these
percolating strips at a particular $V_{N}^{i}$ would be given by
$\Delta x_{i} =W f(V_N^{i}) \Delta V_{N}$. The total conductance
of the specimen system would equal the sum of
the conductances of all these strips. For simplicity, assume that
all the strips have the same mobility. Then, $\sigma_{xx} =
(e\mu/(1+(\mu B)^2)(1/W)(\sum n_{i} \Delta x_{i})$. In
neutrality voltage space, after separating out the
electron and hole components, one obtains $\sigma_{xx} =(e\mu/(1+(\mu
B)^2)[C\int_{\nu-w}^{V_G} (V_{G} - V_{N}) f(V_{N})dV_{N} +
C\int_{V_{G}}^{\nu + w} (V_{N} - V_{G}) f(V_{N})dV_{N}]$.
Similarly, $\sigma_{xy} =(e\mu^{2} B/(1+(\mu
B)^2)[-C\int_{\nu-w}^{V_G} (V_{G} - V_{N}) f(V_{N})dV_{N} +
C\int_{V_{G}}^{\nu + w} (V_{N} - V_{G}) f(V_{N})dV_{N}]$.

Using this model, we evaluated the diagonal resistance $R_{xx} =
\rho_{xx} (L/W)$ and $dR_{xy}/dB = (R_{xy} (0.5 T)/(0.5 T))$. A
least squares data-fit was carried out with $\nu$ and $w$ as the
fitting parameters. The blue trace in Fig. 4(a) illustrates the $f(V_N)$ that served to
obtain the best data fit shown in Fig. 4(b) and Fig. 4(c). This
$f(V_{N})$ peaks at $\nu \approx 72.8V$ and exhibits a width
$w=36.5 V$. The extracted effective electron and hole densities for the best
fit $f(V_N)$ are also shown in Fig. 4(a). Here, the effective hole density
at a particular gate voltage $V_{G} \le \nu+w$ is $n_h^{eff.} (V_{G}) =
C\int_{V_{G}}^{\nu+w} (V_{N}-V_{G})f(V_N) dV_{N}$. Similarly, the effective electron density at  $V_{G} \ge \nu-w$ is $n_e^{eff.}(V_{G}) = C\int_{\nu -
w}^{V_G} (V_{G}-V_{N})f(V_N) dV_{N}$. The notable features are that (a) for a broad range of
$V_G$ about $\nu = 72.8 V$, one finds non-vanishing $n_e^{eff.}$ and
$n_h^{eff.}$, unlike the expectations for $n_{e}$ and $n_{h}$ indicated by the black lines for
ideal graphene, and (b) at $\nu = 72.8 V$, the effective densities are given by $n_h^{eff.} = n_e^{eff.} = 0.49
\times 10^{12} cm^{-2}$. Thus, this approach helps us to determine the effective residual densities of electrons and holes at the nominal neutrality point. Fig. 4(b) shows the experimental $dR_{xy}/dB$ data along with the
best fit, while Fig. 4(c) shows $R_{xx}$ data along with the best
fit. These good fits suggest that the model captures the
transport behavior observed in the measurements.

When the Fermi level in ideal monolayer graphene is brought
towards the Dirac neutrality point, naively, one expects $n_q$ to
progressively vanish, leading to a divergent $dR_{xy}/dB$, see
brown trace Fig. 4(b) (inset), with a discontinuity (dashed red
line in Fig. 4(b), inset) in $dR_{xy}/dB$ vs. $V_{G}$ at the
nominal neutrality point.
Measurements for large area disordered CVD graphene devices shown here illustrate  that the
$dR_{xy}/dB$ changes smoothly from positive to negative values
over a broad range of $V_G$, and $dR_{xy}/dB = 0$ at $V_G=\nu$.
Indeed, in experiment, one might set $V_{G} = \nu$ and obtain
vanishing Hall effect, or close to it, as in Fig. 1, with stable and reproducible results.

This modelling suggests ambipolar transport over a wide range of
$V_{G}$, about $V_G = \nu = 72.8V$, see Fig. 4(a), because there is a concurrent non-vanishing effective electron- and hole- density over a broad span of $V_{G}$. The model
provided here is reminiscent of the uniform-carrier-density two
carrier Drude model. However, this model differs
on the point that instead of two carriers at fixed densities,
there is a gradation in the density of both the electrons and the
holes over the sample. If we set the effective hole density to the hole density, i.e., $n_{h}^{eff.} = n_{h}$, and the effective electron density to the electron density, i.e., $n_{e}^{eff.} = n_{e}$, then this model is formally similar to the two carrier
Drude model. Thus, we can use the
predictions of the two carrier ambipolar Drude model to understand
some observations. For example, in the ambipolar two carrier model
with $\mu_h = \mu_e = \mu$ and $\mu B <<1$, $dR_{xy}/dB = (1/e)
(n_h-n_e)/(n_h+n_e)^2$. Thus, $dR_{xy}/dB \propto n_h-n_e$, and
$dR_{xy}/dB$ should vanish at $V_G=\nu$, where $n_h = n_e$, as observed in experiment, see Fig. 3(b) or
Fig. 4(b). A qualitative picture, which suggests streams rather than puddles, for the $V_{G} = \nu$ condition
is exhibited in Fig. 4(c), inset. At $V_G=\nu$, the applied
current is equally split by symmetry between
counter-propagating holes and electrons. In a magnetic field, the
Hall electric field, $E_H$, within these hole- and electron-
current domains will be oppositely directed, which will lead to a
vanishing global Hall effect.\cite{24,23} Yet, the diagonal
resistance/resistivity would remain finite due to the finite residual densities and the non-vanishing mobility.  For the experimental parameters, this residual
resistivity turns out to be $\approx h/4e^{2}$.\cite{9} Here, the
observed weak temperature dependence of the residual resistivity
is attributed to weak localization and the relative temperature
insensitivity of the transport parameters at the nominal
neutrality point.\cite{30} Fig. 4(c), inset, also helps to qualitatively understand small deviations from
$V_G=\nu$:  Deviation towards the hole regime, i.e, $V_G \le \nu$,
will increase the hole current at the expense of the electron
current. Then, the hole Hall electric field will exceed the
electron Hall electric field, yielding a net hole-like Hall
effect. On the other hand, deviating towards the electron regime,
i.e., $V_{G} \ge \nu$, will increase the electron current at the
expense of the hole current, yielding a net electron-like Hall
effect.

\begin{comment}
Note,
however, that only the fit extracted parameters well away from
$V_{G} = \nu$ should be interpreted within the framework of
ref.\cite{21b}.
\end{comment}
\begin{comment}
In summary, we have experimentally examined the reduced or compensated Hall effect
over a wide range of gate-voltages on either side of the
gate-induced n$\leftrightarrow$p transition in CVD graphene,  along with the diagonal resistance through the Dirac peak which exhibits the $\approx h/4e^2$ residual
resistivity. A distribution of neutrality potentials has been invoked in conjunction with multi-carrier conduction to model the experimental observations. From the modeling, we extracted the effective residual electron- and hole- densities  around the nominal neutrality point. The results indicate co-existing electrons and holes broadly about the nominal neutrality point in large area zero band-gap CVD graphene devices, which suggests domain confined ambipolar currents over the gate-induced n $\leftrightarrow$ p transition.

Curiously, the suggested current
domains are reminiscent of the current domains\cite{24} proposed
for the radiation-induced zero-resistance states in the microwave
irradiated $GaAs/AlGaAs$ 2D electron system.\cite{23,23b}

\section{Acknowledgements}

Low temperature magneto-transport measurements at Georgia State
University are supported by the U.S. Department of Energy, Office
of Basic Energy Sciences, Material Sciences and Engineering
Division under DE-SC0001762. Additional support is provided by the
ARO under W911NF-14-02-0076. Raman measurements were carried out
at the Georgia Institute of Technology in facilities supported by
the MRSEC Program of the National Science Foundation No. DMR
0820382
\pagebreak
\begin{comment}
\section{Methods}

\subsection{CVD Graphene growth}

Commercially available 25$\mu$m thick foils were cut into 2" x 1"
strips, cleaned, and loaded into a 1" quartz tube furnace. The
quartz tube with foil was heated to 1000$^{0}C$  over a period of
1/2 hour in 10 sccm $H_{2}N_{2}$ at a pressure of 250 $\mu mHg$.
The foils were then maintained under these conditions at
$1000^{0}C$ for periods ranging from 15 minutes to 1 hour. At the
end of this annealing period, a flow of $CH_{4}$ at 40 sccm,
supplemented the ongoing $H_{2}N_{2}$ flow. The pressure in the
quartz tube increased to 350 $\mu mHg$ during this graphene growth
phase. The $CH_{4}$ flow was maintained for a period of 1/2 hour
at $1000^{0}C$. Then, the $CH_{4}$ was switched off and the
furnace was rapidly cooled to room temperature.

\subsection{Graphene Transfer}

PMMA was spun onto the top-side of the copper foil
with graphene and baked at 100$^{0}$ C for 1/2 hour. The backside
of the foil was oxygen plasma etched to remove the backside
graphene layer. Then, the foil with PMMA was floated on a ferric
nitrate solution to etch away the copper foil. After foil etch,
the ferric nitrate was replaced with deionized water, and the
deionized water was changed several times. Finally, the
PMMA/Graphene in deionized water was heated to roughly 50$^{0}$ C
for a period of 1/2 hour to complete the cleaning process. After
cool-down, the PMMA with graphene was scooped up onto a 300nm
$SiO_{2}/Si$ wafer and allowed to air dry for a period of 24
hours. After 24 hours, the PMMA was stripped from the $SiO_{2}/Si$
with standard solvents, and blown dry with nitrogen gas.

\subsection{Sample preparation}

Standard  semiconductor photolithography techniques were utilized
to prepare graphene Hall bars on the $SiO_{2}/Si$ chips. The
graphene was mesa-etched using a rf-oxygen plasma. Cr/Au contacts
served to electrically access the CVD graphene Hall bars. The
relevant Hall bar device length (L) between diagonal voltage
contacts to device width (W) ratio is one, and $W = 0.4 mm$.
Back-gating was carried out by biasing the Si substrate with
respect to the graphene layer.

\subsection{Transport measurements}

$SiO_2/Si$ chips with graphene Hall bars were mounted on
chip-carriers. Wires were wedge-bonded to the metallization make
electrical contacts to graphene Hall bars on $SiO_{2}/Si$ chips. A
sample was loaded onto a sample holder and lowered into a liquid
helium cryostat/ superconducting magnet system. Four terminal
electrical measurements were then carried out using standard low
frequency lock-in techniques in a transverse magnetic field.
\end{comment}


\section{Introduction}
Let us consider the quasilinear Neumann system
$$\left\{
\begin{array}{lll}
-\Delta_p u+a(x)|u|^{p-2}u=\lambda c(x) F_u(u,v)& {\rm in} & \Omega,\\
-\Delta_q v+b(x)|v|^{q-2}v=\lambda c(x) F_v(u,v)& {\rm in} & \Omega,\\
\frac{\partial u}{\partial \nu}=\frac{\partial v}{\partial \nu}=0 &
{\rm on} & \partial\Omega,
\end{array}\right.
 \eqno{(N^{p,q}_\lambda)}
 $$
 where $p,q>1;$ $\Omega\subset \mathbb R^N$ is a smooth open bounded domain; $\nu$ denotes the outward unit normal to $\partial \Omega$;
 $a,b,c\in L^\infty(\Omega)$ are some functions;  $\lambda\geq 0$ is a
 parameter; and $F_u$ and $F_v$ denote the partial derivatives of $F\in C^1(\mathbb{R}^2,\mathbb{R})$ with
 respect to the first and second variables, respectively.

Recently, problem $(N^{p,q}_\lambda)$ has been considered by several
authors. For instance, under suitable assumptions on $a,b,c$ and
$F$,  El Manouni and Kbiri Alaoui \cite{EM-KA} proved the existence
of an interval $A\subset (0,\infty)$ such that $(N^{p,q}_\lambda)$
has at least three solutions whenever $\lambda\in A$ and $p,q>N$.
Lisei and Varga \cite{LV} also established the existence of at least
three solutions for the system $(N^{p,q}_\lambda)$ with
nonhomogeneous and nonsmooth Neumann boundary conditions. Di Falco
\cite{DiFalco} proved the existence of infinitely many solutions for
$(N^{p,q}_\lambda)$ when the nonlinear function $F$ has a suitable
oscillatory behavior. Systems similar to $(N^{p,q}_\lambda)$ with
the Dirichlet boundary conditions were also considered by Afrouzi
and Heidarkhani \cite{AH-1, AH-2}, Boccardo and de Figueiredo
\cite{deFig-Boccardo}, Heidarkhani and Tian \cite{HT}, Li and Tang
\cite{LT}, see also references therein.

The aim of the present paper is to describe a new phenomenon for
Neumann systems when the nonlinear term has a subquadratic growth.
In order to avoid technicalities, instead of the quasilinear system
$(N^{p,q}_\lambda)$, we shall consider the semilinear problem
$$\left\{
\begin{array}{lll}
-\Delta u+a(x)u=\lambda c(x) F_u(u,v)& {\rm in} & \Omega,\\
-\Delta v+b(x)v=\lambda c(x) F_v(u,v)& {\rm in} & \Omega,\\
\frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=0 & {\rm
on} & \partial\Omega.
\end{array}\right.
 \eqno{(N_\lambda)}
 $$

We assume that the nonlinear term $F\in
C^1(\mathbb{R}^2,\mathbb{R})$ satisfies the following properties:
\begin{itemize}
  \item[{ $(\bf F_+)$}] $F(s,t)\geq 0$ for every $(s,t)\in \mathbb
  R^2$, $F(0,0)=0,$ and $F\not\equiv 0;$
  \item[{ $(\bf F_0)$}]
  $\lim_{(s,t)\to (0,0)}\frac{F_s(s,t)}{{|s|}+{|t|}}=\lim_{(s,t)\to (0,0)}\frac{F_t(s,t)}{{|s|}+{|t|}}=0;$
  \item[{ $(\bf F_\infty)$}]
   $\lim_{|s|+|t|\to \infty}\frac{F_s(s,t)}{{|s|}+{|t|}}=\lim_{|s|+|t|\to \infty}\frac{F_t(s,t)}{{|s|}+{|t|}}=0.$
\end{itemize}

\begin{example}\rm\label{pelda}
A typical nonlinearity which fulfils hypotheses $(\bf F_+)$, $(\bf
F_0)$ and $(\bf F_\infty)$ is $F(s,t)=\ln(1+s^2t^2).$
\end{example}
We also introduce the set $$\Pi_+(\Omega)=\{a\in L^\infty(\Omega):
{\rm essinf}_\Omega a>0\}.$$ For $a,b,c\in \Pi_+(\Omega)$ and for
$F\in C^1(\mathbb R^2,R)$ which fulfils the hypotheses $(\bf F_+)$,
$(\bf F_0)$ and $(\bf F_\infty)$, we define  the numbers
$$s_F=2\|c\|_{L^1}\max_{(s,t)\neq
(0,0)}\frac{F(s,t)}{\|a\|_{L^1}{s^2}+\|b\|_{L^1}t^2},\ {\rm and}\
S_F=\max_{(s,t)\neq
(0,0)}\frac{|sF_s(s,t)+tF_t(s,t)|}{{\|c/a\|_{L^\infty}^{-1}s^2}+{\|c/b\|_{L^\infty}^{-1}t^2}}.$$
Note that these numbers are finite, positive and $S_F\geq s_F$, see
Proposition \ref{elso-prop} (here and in the sequel,
$\|\cdot\|_{L^p}$ denotes the usual norm of the Lebesgue space
$L^p(\Omega)$, $p\in [1,\infty]$). Our main result reads as follows.

\begin{theorem}\label{fotetel-NA} Let  $F\in C^1(\mathbb{R}^2,\mathbb{R})$ be
a function which satisfies $(\bf F_+)$, $(\bf F_0)$ and $(\bf
F_\infty)$, and $a,b,c\in \Pi_+(\Omega)$. Then, the following
statements hold.
\begin{enumerate}
\item[{\rm (i)}] For every
$0\leq\lambda<S_F^{-1}$, system $(N_\lambda)$ has only the trivial
pair of solution.
\item[{\rm (ii)}] For every $\lambda>s_F^{-1}$, system $(N_\lambda)$ has at least two
distinct, nontrivial pairs of solutions
$(u_\lambda^i,v_\lambda^i)\in H^1(\Omega)^2$, $i\in \{1,2\}$.
 \end{enumerate}
\end{theorem}

\begin{remark}\rm (a) A natural question arises which is still open: how many solutions exist for $(N_\lambda)$  when
$\lambda\in [S_F^{-1},s_F^{-1}]$? Numerical experiments show that
$s_F$ and $S_F$ are usually not far from each other, although their
origins are independent. For instance, if $a=b=c$, and $F$ is from
Example \ref{pelda}, we have $s_F\approx 0.8046$ and $S_F=1.$

(b) Assumptions $(\bf F_+)$, $(\bf F_0)$ and $(\bf F_\infty)$ imply
that there exists $c>0$ such that
\begin{equation}\label{F-re-vonatkozo-szukvadratikus}
    0\leq F(s,t)\leq c(s^2+t^2) \ {\rm for\ all}\ (s,t)\in \mathbb R^2,
\end{equation}
i.e., $F$ has a subquadratic growth. Consequently, Theorem
\ref{fotetel-NA} completes the results of several papers where $F$
fulfils the Ambrosetti-Rabinowitz condition, i.e., there exist
$\theta>2$ and $r>0$ such that
\begin{equation}\label{AR}
0<\theta F(s,t)\leq s F_s(s,t)+t F_t(s,t)\ {\rm for\ all}\
|s|,|t|\geq r.
\end{equation}
Indeed, (\ref{AR}) implies that for some $C_1,C_2>0$, one has
$F(s,t)\geq C_1(|s|^\theta+|t|^\theta)$ for all $|s|,|t|> C_2.$
\end{remark}

The next section contains some auxiliary notions and results, while
in Section \ref{sect-bizonyitas} we prove Theorem \ref{fotetel-NA}.
First, a direct calculation proves (i), while a very recent three
critical points result of  Ricceri \cite{Ri-1} provides the proof of
(ii).

\section{Preliminaries}

A solution for $(N_\lambda)$ is a pair $(u,v)\in H^1(\Omega)^2$ such
that
\begin{equation}\label{va}
    \left\{
\begin{array}{lll}
\int_\Omega (\nabla u \nabla \phi+a(x)u\phi)dx=\lambda \int_\Omega c(x) F_u(u,v)\phi dx& {\rm for\ all} & \phi\in H^1(\Omega),\\
\int_\Omega (\nabla v \nabla \psi+b(x)v\psi)dx=\lambda \int_\Omega
c(x) F_v(u,v)\psi dx& {\rm for\ all} & \psi\in H^1(\Omega).
\end{array}\right.
\end{equation}

Let $a,b,c\in \Pi_+(\Omega).$ We associate to the system
$(N_\lambda)$ the energy functional $I_\lambda:H^1(\Omega)^2\to
\mathbb R$ defined by
$$I_\lambda(u,v)=\frac{1}{2}(\|u\|_a^2+\|v\|_b^2)-\lambda \mathcal F(u,v),$$
where $$\|u\|_a=\left(\int_\Omega |\nabla
u|^2+a(x)u^2\right)^{1/2};\ \|v\|_b=\left(\int_\Omega |\nabla
v|^2+b(x)v^2\right)^{1/2},$$ and $$\mathcal F(u,v)=\int_\Omega
c(x)F(u,v).$$ It is clear that $\|\cdot \|_a$ and $\|\cdot \|_b$ are
equivalent to the usual norm on $H^1(\Omega).$ Note that if $F\in
C^1(\mathbb R^2,R)$  verifies the hypotheses $(\bf F_0)$ and $(\bf
F_\infty)$ (see also relation
(\ref{F-re-vonatkozo-szukvadratikus})), the functional $I_\lambda$
is well-defined, of class $C^1$ on $H^1(\Omega)^2$ and its critical
points are exactly the solutions for $(N_\lambda)$. Since
$F_s(0,0)=F_t(0,0)=0$ from $(\bf F_0)$, $(0,0)$ is a solution of
$(N_\lambda)$ for every $\lambda\geq 0.$

In order to prove Theorem \ref{fotetel-NA} (ii), we must find
critical points for $I_\lambda.$ In order to do this, we recall the
following Ricceri-type three critical point theorem. First, we need
the following notion: if $X$ is a Banach space, we denote by
$\mathcal W_X$ the class of those functionals $E:X\to \mathbb{R}$
that possess the property that if $\{u_n\}$ is a sequence in $X$
converging weakly to $u\in X$ and $\liminf_n E(u_n)\leq E(u)$ then
$\{u_n\}$ has a subsequence strongly converging to $u$.

\begin{theorem} {\rm \cite[Theorem 2]{Ri-1}}\label{bonanno-tetel}
Let $X$ be a separable and reflexive real Banach space, let
$E_1:X\to \mathbb{R}$ be a coercive, sequentially weakly lower
semicontinuous $C^1$ functional belonging to $\mathcal W_X$, bounded
on each bounded subset of $X$ and whose derivative admits a
continuous inverse on $X^*$; and $E_2:X\to \mathbb{R}$ a $C^1$
functional with a compact derivative. Assume that $E_1$ has a strict
local minimum $u_0$ with $E_1(u_0)=E_2(u_0)=0$. Setting the numbers
\begin{equation}\label{tau}
\tau=\max\left\{0,\limsup_{\|u\|\to
\infty}\frac{E_2(u)}{E_1(u)},\limsup_{u\to
u_0}\frac{E_2(u)}{E_1(u)}\right\},
\end{equation}
\begin{equation}\label{chi}
    \chi=\sup_{E_1(u)>0}\frac{E_2(u)}{E_1(u)},
\end{equation}
assume that $\tau<\chi.$

Then, for each compact interval $[a,b]\subset (1/\chi,1/\tau)$
$($with the
 conventions
$1/0=\infty$ and $1/\infty=0$$)$ there exists  $\kappa>0$ with the
following property: for every $\lambda\in [a,b]$ and every $C^1$
functional $E_3:X\to \mathbb{R}$ with a compact derivative, there
exists $\delta>0$ such that for each $\theta\in [0,\delta],$ the
equation $$E_1'(u)-\lambda E_2'(u)-\theta E_3'(u)=0$$ admits at
least three solutions in $X$ having norm less than $\kappa.$
\end{theorem}

\noindent We conclude this section with an  observation which
involves the constants $s_F$ and $S_F$.

\begin{proposition}\label{elso-prop} Let  $F\in C^1(\mathbb{R}^2,\mathbb{R})$ be
a function which satisfies $(\bf F_+)$, $(\bf F_0)$ and $(\bf
F_\infty)$, and $a,b,c\in \Pi_+(\Omega)$. Then the numbers $s_F$ and
$S_F$ are finite, positive and $S_F\geq s_F.$
\end{proposition}

{\it Proof.} It follows by $(\bf F_0)$ and $(\bf F_\infty)$ and by
the continuity of the functions $(s,t)\mapsto
\frac{F_s(s,t)}{|s|+|t|}$, $(s,t)\mapsto \frac{F_t(s,t)}{|s|+|t|}$
away from $(0,0)$,  that there exists $M>0$ such that
$$|F_s(s,t)|\leq M(|s|+|t|)\ {\rm and}\ |F_t(s,t)|\leq M(|s|+|t|)\ {\rm for\ all}\ (s,t)\in \mathbb
R^2.$$ Consequently, a standard mean value theorem together with
$(\bf F_+)$ implies that
\begin{equation}\label{M-es-becsles}
    0\leq F(s,t)\leq 2M(s^2+t^2)\ {\rm for\ all}\ (s,t)\in \mathbb R^2.
\end{equation}
 We now prove that
\begin{equation}\label{zeroban-F}
    \lim_{(s,t)\to (0,0)}\frac{F(s,t)}{{s^2}+{t^2}}=0\ {\rm and}\ \lim_{|s|+|t|\to \infty}\frac{F(s,t)}{{s^2}+{t^2}}=0.
\end{equation}

By $(\bf F_0)$ and $(\bf F_\infty)$, for every $\eps>0$ there exists
$\delta_\eps\in (0,1)$ such that for every $(s,t)\in \mathbb R^2$
with $|s|+|t|\in (0,\delta_\eps)\cup (\delta_\eps^{-1},\infty)$, one
has
\begin{equation}\label{eps-F}
    \frac{|F_s(s,t)|}{{|s|}+{|t|}}<\frac{\eps}{4}\ {\rm and}\
\frac{|F_t(s,t)|}{{|s|}+{|t|}}<\frac{\eps}{4}.
\end{equation}
By (\ref{eps-F}) and the mean value theorem, for every $(s,t)\in
\mathbb R^2$ with $|s|+|t|\in (0,\delta_\eps)$,
 we have
\begin{eqnarray*}
  F(s,t) &=& F(s,t)-F(0,t)+F(0,t)-F(0,0) \\
   &\leq& \frac{\eps}{2}(s^2+t^2)
\end{eqnarray*}
which gives the first limit in (\ref{zeroban-F}). Now, for every
$(s,t)\in \mathbb R^2$ with
$|s|+|t|>\delta_\eps^{-1}\max\{1,\sqrt{8M/\eps}\}$, by using
(\ref{M-es-becsles}) and (\ref{eps-F}),
 we have
\begin{eqnarray*}
  F(s,t) &=& F(s,t)-
  F\left(\frac{\delta_\eps^{-1}}{|s|+|t|}s,t\right)
  +F\left(\frac{\delta_\eps^{-1}}{|s|+|t|}s,t\right)
  -F\left(\frac{\delta_\eps^{-1}}{|s|+|t|}s,\frac{\delta_\eps^{-1}}{|s|+|t|}t\right)\\ &&+F\left(\frac{\delta_\eps^{-1}}{|s|+|t|}s,\frac{\delta_\eps^{-1}}{|s|+|t|}t\right)  \\
   &\leq& \frac{\eps}{4}(|s|+|t|)^2+2M \delta_\eps^{-2}\\&\leq&
   \eps(s^2+t^2),
\end{eqnarray*}
which leads us to the second limit in (\ref{zeroban-F}).

The facts above show that the numbers $s_F$ and $S_F$ are finite.
Moreover, $s_F>0$. We now prove that $S_F\geq s_F.$ To do this, let
$(s_0,t_0)\in \mathbb R^2\setminus \{(0,0)\}$ be a maximum point of
the function
$(s,t)\mapsto\frac{F(s,t)}{\|a\|_{L^1}{s^2}+\|b\|_{L^1}t^2}.$ In
particular, its partial derivatives vanishes at $(s_0,t_0)$,
yielding
$$F_s(s_0,t_0)(\|a\|_{L^1}{s_0^2}+\|b\|_{L^1}t_0^2)=2\|a\|_{L^1}s_0F(s_0,t_0);$$
$$F_t(s_0,t_0)(\|a\|_{L^1}{s_0^2}+\|b\|_{L^1}t_0^2)=2\|b\|_{L^1}t_0F(s_0,t_0).$$
From the two relations above we obtain that
$$s_0 F_s(s_0,t_0)+t_0 F_t(s_0,t_0)=2F(s_0,t_0).$$
On the other hand, since $a,b,c\in \Pi_+(\Omega)$, we have that
$$\|c\|_{L^1}=\int_\Omega c(x)dx=\int_\Omega \frac{c(x)}{a(x)}
a(x)dx\leq \left\|\frac{c}{a}\right\|_{L^\infty}\int_\Omega
a(x)dx=\left\|\frac{c}{a}\right\|_{L^\infty}\|a\|_{L^1},$$ thus
$\|c/a\|_{L^\infty}^{-1}\leq \|a\|_{L^1}/\|c\|_{L^1}$  and in a
similar way $\|c/b\|_{L^\infty}^{-1}\leq \|b\|_{L^1}/\|c\|_{L^1}$.
Combining these inequalities with the above argument,
we conclude that $S_F\geq s_F.$  \hfill $\square$\\

\vspace{-0.5cm}
\section{Proof of Theorem \ref{fotetel-NA}}\label{sect-bizonyitas}

In this section we assume that the assumptions of Theorem
\ref{fotetel-NA} are fulfilled.\\

\noindent {\bf Proof of Theorem \ref{fotetel-NA} (i).} Let $(u,v)\in
H^1(\Omega)^2$ be a solution of $(N_\lambda).$ Choosing $\phi=u$ and
$\psi=v$ in (\ref{va}),  we obtain that
\begin{eqnarray*}
  \|u\|_a^2 + \|v\|_b^2&=& \int_\Omega (|\nabla
u|^2+a(x)u^2+|\nabla
v|^2+b(x)v^2) \\
   &=& \lambda \int_\Omega c(x)(F_u(u,v)u+F_v(u,v)v)
 \\
   &\leq & \lambda S_F \int_\Omega c(x)({\|c/a\|_{L^\infty}^{-1}u^2}+{\|c/b\|_{L^\infty}^{-1}v^2}) \\
   &\leq &\lambda S_F \int_\Omega (a(x)u^2+b(x)v^2)\\ &\leq &\lambda S_F
   (\|u\|_a^2 + \|v\|_b^2).
\end{eqnarray*}
Now, if $0\leq \lambda< S_F^{-1},$ we necessarily have that
$(u,v)=(0,0)$, which concludes the proof.
\newpage

\noindent {\bf Proof of Theorem \ref{fotetel-NA} (ii).} In Theorem
\ref{bonanno-tetel} we choose $X=H^1(\Omega)^2$ endowed with the
norm $\|(u,v)\|=\sqrt{\|u\|_a^2+\|v\|_b^2}$, and $E_1,
E_2:H^1(\Omega)^2\to \mathbb R$ defined by
$$E_1(u,v)=\frac{1}{2}\|(u,v)\|^2\ {\rm and}\ E_2(u,v)=\mathcal F(u,v).$$
It is clear that both $E_1$ and $E_2$ are $C^1$ functionals and
$I_\lambda=E_1-\lambda E_2.$ It is also a standard fact that $E_1$
is a coercive, sequentially weakly lower semicontinuous functional
which belongs to $\mathcal W_{H^1(\Omega)^2}$, bounded on each
bounded subset of $H^1(\Omega)^2$, and its derivative admits a
continuous inverse on $(H^1(\Omega)^2)^*.$ Moreover, $E_2$ has a
compact derivative since $H^1(\Omega)\hookrightarrow L^p(\Omega)$ is
a compact embedding for every $p\in (2,2^*)$.

Now, we prove that the functional $(u,v)\mapsto
\frac{E_2(u,v)}{E_1(u,v)}$ has similar properties as the function
$(s,t)\mapsto \frac{F(s,t)}{s^2+t^2}.$ More precisely, we shall
prove that
\begin{equation}\label{egy-ketto}
    \lim_{\|(u,v)\|\to 0}\frac{E_2(u)}{E_1(u)}=\lim_{\|(u,v)\|\to \infty}\frac{E_2(u)}{E_1(u)}=
    0.
\end{equation}
First, relation
 (\ref{zeroban-F}) implies that for every $\eps>0$ there exists
$\delta_\eps\in (0,1)$ such that for every $(s,t)\in \mathbb R^2$
with $|s|+|t|\in (0,\delta_\eps)\cup (\delta_\eps^{-1},\infty)$, one
has
\begin{equation}\label{eps-F-uj}
    0\leq \frac{F(s,t)}{{s^2}+t^2}<\frac{\eps}{4\max\{\|c/a\|_{L^\infty},\|c/b\|_{L^\infty}\}}.
\end{equation}
Fix  $p\in (2,2^*).$ Note that the continuous function $(s,t)\mapsto
\frac{F(s,t)}{|s|^p+|t|^p}$ is bounded on the set $\{(s,t)\in
\mathbb R^2:|s|+|t|\in [\delta_\eps,\delta_\eps^{-1}]\}$. Therefore,
for some $m_\eps>0,$ we have that in particular
$$0\leq F(s,t)\leq \frac{\eps}{4\max\{\|c/a\|_{L^\infty},\|c/b\|_{L^\infty}\}}(s^2+t^2) + m_\eps (|s|^p+|t|^p)\ \ {\rm for\ all}\ (s,t)\in \mathbb R^2.$$
Therefore, for each $(u,v)\in H^1(\Omega)^2,$ we get
\begin{eqnarray*}
 0\leq E_2(u,v) &=& \int_\Omega c(x)F(u,v)  \\
   &\leq& \int_\Omega c(x)\left[\frac{\eps}{4\max\{\|c/a\|_{L^\infty},\|c/b\|_{L^\infty}\}}(u^2+v^2) + {m_\eps} (|u|^p+|v|^p)\right] \\
   &\leq& \int_\Omega \left[\frac{\eps}{4}(a(x)u^2+b(x)v^2) + {m_\eps} c(x)(|u|^p+|v|^p)\right] \\
   &\leq& \frac{\eps}{4}\|(u,v)\|^2+ {m_\eps} \|c\|_{L^\infty}S_{p}^{p}(\|u\|_a^{p}+\|v\|_b^{p})\\
   &\leq& \frac{\eps}{4}\|(u,v)\|^2+ {m_\eps} \|c\|_{L^\infty}S_{p}^{p}\|(u,v)\|^p,
\end{eqnarray*}
where $S_{l}>0$ is the best constant in the inequality
$\|u\|_{L^l}\leq S_l \min \{\|u\|_a,\|u\|_b\}$ for every $u\in
H^1(\Omega)$, $l\in (1,2^*)$ (we used the fact that the function
$\alpha\mapsto (s^\alpha+t^\alpha)^\frac{1}{\alpha}$ is decreasing,
$s,t\geq 0$). Consequently, for every $(u,v)\neq (0,0),$ we obtain
\begin{eqnarray*}
  0\leq \frac{E_2(u,v)}{E_1(u,v)} &\leq& \frac{\eps}{2}+ {2m_\eps}
  \|c\|_{L^\infty}S_{p}^{p}\|(u,v)\|^{p-2}.
\end{eqnarray*}
Since $p>2$ and $\eps>0$ is arbitrarily small when $(u,v)\to 0,$ we
obtain the first limit from (\ref{egy-ketto}).

 Now, we fix $r\in (1,2).$ The continuous function $(s,t)\mapsto
\frac{F(s,t)}{|s|^r+|t|^r}$ is bounded on the set $\{(s,t)\in
\mathbb R^2:|s|+|t|\in [\delta_\eps,\delta_\eps^{-1}]\}$, where
$\delta_\eps\in (0,1)$ is from (\ref{eps-F-uj}). Combining this fact
with (\ref{eps-F-uj}), one can find a number $M_\eps>0$ such that
$$0\leq F(s,t)\leq \frac{\eps}{4\max\{\|c/a\|_{L^\infty},\|c/b\|_{L^\infty}\}}(s^2+t^2) + M_\eps (|s|^r+|t|^r)\ \ {\rm for\ all}\ (s,t)\in \mathbb R^2.$$
The H\"older inequality and a similar calculation as above show that
\begin{eqnarray*}
  0\leq E_2(u,v) &\leq &\frac{\eps}{4}\|(u,v)\|^2+ 2^{1-\frac{r}{2}}{M_\eps} \|c\|_{L^\infty}S_{r}^{r}\|(u,v)\|^{r}.
\end{eqnarray*}
For every $(u,v)\neq (0,0)$, we have that
\begin{eqnarray*}
  0\leq \frac{E_2(u,v)}{E_1(u,v)} &\leq& \frac{\eps}{2}+ {2^{2-\frac{r}{2}}M_\eps}
  \|c\|_{L^\infty}S_{r}^{r}\|(u,v)\|^{r-2}.
\end{eqnarray*}
 Due to the arbitrariness of $\eps>0$ and $r\in (1,2)$, by letting the limit $\|(u,v)\|\to
 \infty$,  we obtain the second relation from (\ref{egy-ketto}).

Note that $E_1$ has a strict global minimum $(u_0,v_0)=(0,0)$, and
$E_1(0,0)=E_2(0,0)=0.$ The definition of the number $\tau$ in
Theorem \ref{bonanno-tetel}, see (\ref{tau}), and the limits in
(\ref{egy-ketto}) imply that $\tau=0.$ Furthermore, since
$H^1(\Omega)$ contains the constant functions on $\Omega$, keeping
the notation from (\ref{chi}), we obtain
$$\chi=\sup_{E_1(u,v)>0}\frac{E_2(u,v)}{E_1(u,v)}\geq 2\|c\|_{L^1}\max_{(s,t)\neq
(0,0)}\frac{F(s,t)}{\|a\|_{L^1}{s^2}+\|b\|_{L^1}t^2}=s_F.$$
Therefore, applying Theorem \ref{bonanno-tetel} (with $E_3\equiv0$),
we obtain that in particular for every $\lambda\in
(s_F^{-1},\infty)$, the equation $I'_\lambda(u,v)\equiv
E'_1(u,v)-\lambda E_2'(u,v)=0$ admits at least three distinct pairs
of solutions in $H^1(\Omega)^2.$  Due to condition $(\bf F_0)$,
system $(N_\lambda)$ has the solution $(0,0)$. Therefore, for every
$\lambda>s_F^{-1}$, the system $(N_\lambda)$ has at least two
distinct, nontrivial pairs of solutions, which concludes the proof.

\begin{remark}\rm\label{remark-utolso}
The conclusion of Theorem \ref{bonanno-tetel} gives a much more
precise information about the Neumann system $(N_\lambda)$; namely,
one can see that $(N_\lambda)$ is stable with respect to small
perturbations. To be more precise, let us consider the perturbed
system
$$\left\{
\begin{array}{lll}
-\Delta u+a(x)u=\lambda c(x) F_u(u,v)+\mu d(x) G_u(u,v)& {\rm in} & \Omega,\\
-\Delta v+b(x)v=\lambda c(x) F_v(u,v)+\mu d(x) G_v(u,v)& {\rm in} & \Omega,\\
\frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=0 & {\rm
on} & \partial\Omega.
\end{array}\right.
 \eqno{(N_{\lambda,\mu})}
 $$
 where $\mu\in \mathbb R,$ $d\in L^\infty(\Omega)$, and $G\in C^1(\mathbb R^2,\mathbb R)$ is a function such that for some $c>0$ and $1<p<2^*-1,$
 $$\max\{|G_s(s,t)|,|G_t(s,t)|\}\leq c(1+|s|^p+|t|^p)\ {\rm for\ all}\ (s,t)\in \mathbb R^2.$$
 One can prove in a standard manner that $E_3:H^1(\Omega)^2\to \mathbb
 R$ defined by $$E_3(u,v)=\int_{\Omega} d(x)G(u,v)dx,$$
 is of class $C^1$ and it has a compact derivative. Thus, we may apply Theorem \ref{bonanno-tetel} in its generality to show  that
  for small enough values of $\mu$ system $(N_{\lambda,\mu})$ still
 has three distinct pairs of solutions.
\end{remark}