\section{Introduction}

In heavy ion collisions at very high energies,
soft (small-x) partons, mostly
gluons, constitute a big portion of
the wave function of the colliding
nuclei. Since the occupation number of
the gluons is very high, they
can be treated in the framework of
classical Yang-Mills theory. When
two nuclei pass through each other,
the classical nonabelian fields from each
nucleus interact with each other and form color flux tubes between
the receding nuclei, which ultimately lead to the quark gluon plasma.
The energy density distribution of the initial state
determines the multiplicity and spectrum
of the final particles seen in the detectors.
The initial energy density distribution
is an input to the hydrodynamics description
of the evolving quark gluon
plasma and, in principle, it can be calculated
ab initio from Quantum Chromodynamics
(QCD).

The energy initially deposited in the color flux tubes can be related to
the two-point vector potential correlation function
 $\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{b}(\boldsymbol{x}^{\,\prime})\rangle$ where
the setup may be two (transverse) or three dimensional,
depending on whether or not the Lorentz contracted nucleus is assumed to have a longitudinal
thickness in the lab frame. This correlator can be derived analytically from the color
charge density correlator
$\langle\rho^{a}(\boldsymbol{x})\rho^{b}(\boldsymbol{x}^{\,\prime})\rangle$
by using the classical Yang-Mills equations.

The color charge density of a nucleus $\rho^{a}(\boldsymbol{x})$
during a collision cannot be known;
on the other hand, the fluctuations in the
color charge density can be studied in the effective field theory approach
with ensemble averaging.
In the Color Glass Condensate (CGC), a framework for
slowly evolving high density gluons within the ultrarelativistic nucleus,
the fast partons are seen as sources
of the small-x (soft) gluon radiation, where x
is the
fraction of the total longitudinal momentum carried by a parton.
After integrating out the fast partons,
the observables can be calculated by averaging them over the ensemble of all
possible color charge configurations. Once the correlator
$\langle\rho^{a}(\boldsymbol{x})\rho^{b}(\boldsymbol{x}^{\,\prime})\rangle$
is specified,
it can be linked to the vector field correlator
$\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{b}(\boldsymbol{x}^{\,\prime})\rangle$
from which the gluon distribution function $xg(x,Q^{2})$ and other observables,
such as the initial energy density,
can be calculated. Here $\langle \cdots \rangle$
denotes the ensemble average. See Ref. \cite{Iancu:2003xm}
for a comprehensive review.

In the early formulation of the CGC by McLerran and Venugopalan
\cite{McLerran:1993ni,McLerran:1993ka,
McLerran:1994vd,Ayala:1995kg,Ayala:1995hx},
the spectrum of the Gaussian color fluctuations was taken to be
white noise. This results in arbitrarily long wavelength fluctuations where
$\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{b}(\boldsymbol{x}^{\,\prime})\rangle$
diverges in the infrared. This problem originates from the fact that
confinement effects (color neutrality) have not been taken into
account. Adding a gluon mass to the gluon
propagator so as to bypass this problem breaks gauge invariance
which is required later to convert
the solution of the classical Yang-Mills equations
from the axial gauge to the light-cone gauge by means of Wilson lines.

In the original McLerran-Venugopalan (MV) model, the colliding nuclei
are considered
to be two-dimensional infinitely thin sheets traveling at the speed of light.
In this approximation,
the model does not depend on the longitidunal coordinate,
hence there is no dependence on the momentum fraction x
of a given parton.
In other words, all of the partons in the nucleus have
the same momentum fraction x.
Being an artifact of the infinite momentum
frame, the lack of x dependence in the model does not reflect the
true nature of the x dependent gluon distribution functions.

In this paper, we consider a color neutral three dimensional
McLerran-Venugopalan (3dMVn) model
by Lam and Mahlon \cite{Lam:1999wu,Lam:2000nz},
where the spectrum of the Gaussian fluctuations is taken to be infrared-safe
colored noise. Colored noise creates a correlation hole
in the correlation function on the size
of a nucleon and this leads to color neutrality.
(The term ``colored noise'' is not related to the color
charge of QCD.)
 This model comes with another important
 feature that the incoming nuclei are
not exactly on the light cone. Longitudinal coordinate dependence
produces an x dependent gluon distribution function. This enables us to
compare the model with the data parametrization of gluon distribution
functions and therefore fix
the parameters of the 3dMVn model. Once this is achieved,
the 3dMVn model can be used to calculate the initial energy
density distribution in heavy ion collisions.

The MV model provides a framework for calculating
the vector field correlation function
$\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{b}(\boldsymbol{x}^{\,\prime})\rangle$
from which the gluon distribution function as well as the
initial energy density distribution
can be calculated.
We emphasize that our ultimate goal is calculating the latter,
which will be pursued in a follow-up paper.
In this paper we focus on the calculation of gluon distribution function
and its comparison with data for the purpose of fixing
the free parameters
in the 3dMVn model.

The outline of the paper is as follows.
First, we give a brief overview of the color neutral and x-dependent version of the McLerran-Venugopalan model (3dMVn) given by
Lam and Mahlon \cite{Lam:1999wu,Lam:2000nz}. Next, we calculate
the gluon distribution function from the 3dMVn model.
This model works at low-$Q^{2}$
where no data for the gluon distribution functions is available.
For that reason, we evolve our results
with the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation to higher momenta where data is available.
This is followed by a discussion of the data
parametrization that we use to compare with the 3dMVn model.
Lastly, we compare the 3dMVn
with the data parametrization.
The final section summarizes our analysis and includes a
discussion regarding how the results will be used in
the future to calculate initial
energy density distribution.

\section{Color Neutrality}

In the original two dimensional MV model, the average and the fluctuations
of the color charge density of a nucleus are given by
\begin{eqnarray}
\langle\rho^{a}(\boldsymbol{x}_\perp^{~} )\rangle & = & 0,\label{eq:vev}\\
\langle\rho^{a}(\boldsymbol{x}_\perp^{~} )
\rho^{b}(\boldsymbol{x}^{\,\prime}_{\perp})\rangle & = &
\delta^{ab}\mu_{A}^{2} \, {\cal D}({\boldsymbol{x}_\perp^{~} }
-\boldsymbol{x}^{\,\prime}_{\perp}),\label{eq:variance}
\end{eqnarray}
where $\boldsymbol{x}_\perp^{~} $ is the
transverse coordinate system for the
infinitely thin nucleus, $\mu_{A}^{2}$ is the
average color charge density squared
per unit (transverse) area for a nucleus, and
$\langle\cdots\rangle$ denotes the ensemble average.
The function
${\cal D}$ determines the spectrum of the fluctuations.
When white noise
is assumed, ${ {\cal D}(\boldsymbol{x}_\perp^{~}
-\boldsymbol{x}^{\,\prime}_{\perp})=\delta^{2}(\boldsymbol{x}_\perp^{~}
-\boldsymbol{x}^{\,\prime}_{\perp}) }$,
no correlation occurs between different points in a nucleus. This
also means that fluctuations in each momentum mode, including the
zero mode, are equally likely and there is no characteristic scale.
In this case, even though the correlator
$\langle A_{i}^{a}(\boldsymbol{x}_\perp^{~})A_{i}^{b}
(\boldsymbol{x}^{\,\prime}_{\perp})\rangle$
is calculated with an infrared cutoff  $\Lambda_{\rm{QCD}}$
\cite{JalilianMarian:1996xn}, it still diverges like
$(\boldsymbol{x}_\perp^2)^{\boldsymbol{x}_\perp^2}$ in the infrared.
Therefore, the Fourier transform
of it, which
is necessary to calculate the gluon distribution function,
does not exist \cite{Lam:1999wu}.
This cut-off can also be understood as a gluon mass in
the gluon propagator
in the form
$  ({\boldsymbol{q}_\perp^{2}+m_{\rm{gluon}}^{2}})^{-1} $.
It explicitly breaks
the gauge invariance that is needed to convert the solutions of
classical Yang-Mills equations from the axial gauge to
the light-cone gauge
where Wilson lines are to be used.

In reality, a nucleus is color neutral
on scales larger than a
nucleon size. A color neutral correlation function
can be derived from a simple model of a nucleus
where nucleons are composed of quark and antiquark
pairs \cite{Kovchegov:1996ty}.
If we consider a two dimensional nucleus for a moment,
the assumption that only the quark and antiquark pair from the same
nucleon can interact with each other produces a correlator of the
form \cite{Lam:1999wu}
\begin{equation}
\langle\rho^{a}(\boldsymbol{x}_\perp^{~} )\rho^{b}
(\boldsymbol{x}^{\,\prime}_{\perp})\rangle=\delta^{ab}\mu_{A}^{2}
\left[\delta^{2}(\boldsymbol{x}_\perp^{~} -\boldsymbol{x}^{\,\prime}_{\perp})-
\frac{\exp\left[-\vert \boldsymbol{x}_{\perp}^{~}
-\boldsymbol{x}^{\,\prime}_{\perp}\vert^{2}
/ 4 \lambda^{2} \right]}{4\pi \lambda^{2}}
\right].\label{eq:colorednoise}
\end{equation}
Here the parameter $\lambda$, which
we will refer to as a correlation length,
is of the order of
a nucleon size $\sim1\,\rm{fm}$.
This parameter takes away the need for a sharp
infrared cutoff as it is used in the original MV treatment; moreover,
it makes the Fourier transform of the vector field correlation function
well-defined and it renders the model infrared safe.
In Eq. (\ref{eq:colorednoise}), decorrelation
due to the white noise ${\cal D}(\boldsymbol{x}_\perp^{~}
-\boldsymbol{x}^{\,\prime}_{\perp})=\delta^{2}(
\boldsymbol{x}_\perp^{~} -\boldsymbol{x}^{\,\prime}_{\perp})$
is modified by the last term for distances
$\vert \boldsymbol{x}_\perp^{~} -\boldsymbol{x}^{\,\prime}_{\perp}\vert
\lesssim \lambda$,
reflecting the assumption that there is a correlation between
the partons confined to a region smaller than the nucleon size.
The second part in Eq. (\ref{eq:colorednoise}) removes the zero mode
$\vert \boldsymbol{q}_{\perp}^{~} \vert=0$ from
the white noise spectrum, hence colored noise.
The Fourier transform of the
correlator in Eq. (\ref{eq:colorednoise}) is
proportional to
$1-e^{-\lambda^{2}\boldsymbol{q}_{\perp}^{2}}$,
which vanishes as
$\vert \boldsymbol{q}_{\perp}^{~} \vert \rightarrow0$.
This makes the model infrared safe.

In this paper, we adopt a three dimensional correlator
given by Lam and Mahlon \cite{Lam:2000nz}
\begin{equation}
\langle\rho^{a}(0,0)\rho^{b}(\boldsymbol{x})\rangle=\delta^{ab}
\kappa_{A}^{3}
\left[\delta^{3}(\boldsymbol{x})-\frac{3}{4\pi \lambda^{2}}
\frac{
\exp \left( - \sqrt{3} \vert \boldsymbol{x} \vert / \lambda \right)
}{
\vert \boldsymbol{x} \vert
}
\right].\label{eq:Yukawa}
\end{equation}
where $\kappa_{A}^{3}$ is the three dimensional
average color charge density
\begin{equation}
\kappa_{A}^{3}=\frac{3 A C_F}{N_{c}^{2}-1}
\frac{1}{V}=\frac{3 A}{2 N_c} \frac{1}{V},
\label{kappa}
\end{equation}
Here $N_c$ is the number of colors, $V$ is the
volume of a nucleus and the color factor
is defined as $C_F=(N_{c}^{2}-1)/(2 N_c)$.
We define $d^{3}\boldsymbol{x}\equiv dx_{\parallel}d^{2}
\boldsymbol{x}_\perp^{~} $
where $x_{\parallel}$ is the longitudinal coordinate in the direction
of the beam axis, and $\boldsymbol{x}_\perp^{~}$
is the coordinate on the transverse
plane. The Fourier transform of the corralator
(\ref{eq:Yukawa})
is given by
\begin{equation}
\widetilde{\langle\rho^{a}\rho^{b}\rangle}=\delta^{ab}
\kappa_{A}^{3}
\left[1-\frac{1}{1+ \lambda^{2} \boldsymbol{q}^{2}/3} \right].
\label{eq:Yukawa-Fouriered}
\end{equation}
Despite the slight difference between the correlators
in Eqs. (\ref{eq:colorednoise})
and (\ref{eq:Yukawa}), they produce
similar results. It can be readily seen from
Eq. (\ref{eq:Yukawa-Fouriered})
that the zero mode does not exist in the spectrum
since $\widetilde{\langle\rho^{a}\rho^{b}\rangle}\rightarrow0$
as $\vert \boldsymbol{q} \vert \rightarrow0$.

\section{Longitudinal Dependence}

In ultrarelativistic heavy ion collisions, the infinite momentum frame
provides a good starting point for the parton picture of the nucleus.
However, the nucleus becomes infinitely thin when it is exactly
on the light cone and the gluon distribution function
turns out to be independent of x. Here we review
the formulation by
Lam and Mahlon \cite{Lam:2000nz} for the
case where the nucleus is boosted to speed $\beta$.
In the lab frame, the thickness of a nucleus of radius $R$
becomes of the order of $R/\gamma$.

The current for
a nucleus moving in the $+z$
direction
is
\begin{equation}
J_{r}^{0}=\rho(-z,{\boldsymbol{x}_\perp^{~} }_{r});
\,\,\,\,\boldsymbol{J}_{r}=0, \label{eq:rest-current1}
\end{equation}
where the subscript ``r'' stands for the rest frame
and the negative sign in front of the
$z$ is for later convenience. With the redefinition
$x^{\pm}=-x_{\mp}=(t\pm z)/\sqrt{2}$,
the current can be rewritten in the light-cone coordinates (still
in the rest frame),
\begin{equation}
J_{r}^{+}=J_{r}^{-}=\frac{1}{\sqrt{2}}\rho\left(
\frac{1}{\sqrt{2}}(x_{r}^{-}-x_{r}^{+}),
{\boldsymbol{x}_\perp^{~} }_{r}\right);\,\,\,\,
{\boldsymbol{J}_\perp}_{r}=0.\label{eq:rest-light-cone}
\end{equation}

When we go from the rest frame to the lab frame where the nucleus
moves with speed $\beta$, the current in Eq. (\ref{eq:rest-light-cone})
becomes
\begin{equation}
J^{+}=\frac{1}{\varepsilon}\rho\left(\frac{1}{\varepsilon}
x^{-}-\frac{\varepsilon}{2}x^{+},
\boldsymbol{x}_\perp^{~}  \right);\,\,\,\,
J^{-}=\frac{\varepsilon}{2}J^{+};\,\,\,\,
{\boldsymbol{J}_\perp}=0,\label{eq:lab-light-cone}
\end{equation}
where
\begin{equation}
\varepsilon=\sqrt{\frac{2(1-\beta)}{1+\beta}}.
\end{equation}
Here we can define a new longitudinal coordinate
\begin{equation}
 x_{\parallel}\equiv\frac{1}{\varepsilon}x^{-}
 -\frac{\varepsilon}{2}x^{+},
\label{newcoord}
\end{equation}
which is essentially the Lorentz
tranformation of $-z=\left( x^{-}-x^{+} \right) / \sqrt{2}$
(see the definition below Eq. (\ref{eq:rest-current1})).
The Eq. (\ref{eq:lab-light-cone})
can be contrasted with the commonly used current
in the infinite momentum
frame (imf) where the nucleus is infinitely thin,
\begin{equation}
J_{\rm{imf}}^{+}=\delta(x^{-})\rho(\boldsymbol{x}_\perp^{~} );\,\,\,\,
J_{\rm{imf}}^{-}=0;\,\,\,\,
{\boldsymbol{J}_\perp}_{\rm{imf}}=0.
\end{equation}

\section{Calculation and Evolution of The Gluon Distribution Function}

We now turn to the calculation of the gluon
distribution function from the 3dMVn model.
Later we will compare it with data
in order to fix the parameters $\alpha_s$ and $\lambda$.
The gluon distribution function of a nucleus with a
baryon number $A$
can be expressed as an integral of the gluon number
density over the transverse
momenta,
\begin{equation}
g_{A}(x,Q^{2})\equiv\int_{\vert  {\boldsymbol{q}_\perp} \vert
\leqslant Q}d^{2}{\boldsymbol{q}_\perp} \frac{dN}{dxd^{2}
\boldsymbol{q}_\perp}.\label{eq:gluon-distro}
\end{equation}
The gluon number density in the lab frame for a nucleus moving with speed
$ $$\beta$ is given as a Fourier transform of the two-point vector
field correlation function \cite{Lam:2000nz}
\begin{equation}
\frac{dN}{dq_{\parallel}d^{2} {\boldsymbol{q}_\perp}  }\equiv
\frac{q_{\parallel}}{4\pi^{3}}\int d^{3}\boldsymbol{x}
\int d^{3}\boldsymbol{x}^{\,\prime}
e^{i \boldsymbol{q} \cdot(\boldsymbol{x}-\boldsymbol{x}^{\,\prime})}
\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{a}(\boldsymbol{x}^{\,\prime})
\rangle,\label{eq:gluon-number-density}
\end{equation}
Here $q_{\parallel}$ is the momentum conjugate to the coordinate
defined in Eq. (\ref{newcoord}).
The Eq. (\ref{eq:gluon-number-density}) can be related to
 $ dN/dxd^{2} \boldsymbol{q}_\perp$
by using the relation $x\equiv q_{\parallel} / m$,
which itself can be derived from the definition
$x\equiv q^{+}/Q^{+}= \varepsilon q^{+}/m$.
Here $m$ is the nucleon mass, $q^{+}$ and $Q^{+}$ are the momenta
of the gluon and the nucleon.

Solving the classical Yang-Mills
equations for the source in Eq. (\ref{eq:rest-light-cone})
in the covariant gauge and transforming the solution into the light-cone
gauge by using Wilson lines, one can express the
correlator $ $$\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{a}
(\boldsymbol{x}^{\,\prime})\rangle$
in Eq. (\ref{eq:gluon-number-density}) in terms of
the color charge density
correlator $ $$\langle\rho^{a}(\boldsymbol{x})\rho^{a}
(\boldsymbol{x}^{\,\prime})\rangle$.
We skip the details of this calculation here and refer the reader
to the original paper \cite{Lam:2000nz}. Using the
correlator in Eq. (\ref{eq:Yukawa}),
the gluon number density can be written as
\begin{equation}
\frac{dN}{dxd^{2} {\boldsymbol{q}_\perp}}=\frac{
8 A
\alpha_{s}}{\pi^{2}}\frac{1}{x}\int d^{2} {\boldsymbol{\Delta}_\perp}
e^{i {\boldsymbol{q}_\perp}\cdot
{\boldsymbol{\Delta}_\perp}}{\cal L}(x;{\boldsymbol{\Delta}_\perp})
{\cal E}(v^{2}L({\boldsymbol{\Delta}_\perp})),\label{eq:gluon-distro-LM}
\end{equation}
Here the integration is over the transverse coordinate ${\boldsymbol{\Delta}_\perp}=
\boldsymbol{x}_\perp^{~}  -\boldsymbol{x}^{\,\prime}_{\perp}$.
The functions $\cal{L}$ and $L$ are convolutions
of the two gluon propagators
at two different points in the same nucleus
connected by the three dimensional noise term
${\cal D}( \boldsymbol{x} - \boldsymbol{x}^{\,\prime})$.
These two functions
can be seen as a pair distribution function.
The nuclear correction factor
$\cal{E}$ takes into account the nuclear
geometry and $v^2$ controls the strength
of dependence on this geometry.

The functions used in Eq. (\ref{eq:gluon-distro-LM})
are given as \cite{Lam:2000nz}
\begin{equation}
{\cal L} (x;{\boldsymbol{\Delta}_\perp})=-\frac{1}{12 \pi}
(x m \lambda)^2 K_0(x m \Delta_\perp),
\end{equation}
and
\begin{equation}
L({\boldsymbol{\Delta}_\perp})=-\frac{\lambda^2}{6 \pi}
\left[ K_0\left(
\frac{\sqrt{3} \Delta_\perp}{\lambda}\right)+
\rm{ln}\left(\frac{\sqrt{3}
\Delta_\perp}{2\lambda}\right)+\gamma_E    \right],
\end{equation}
where $\Delta_\perp=\vert \boldsymbol{\Delta}_\perp \vert$. In the calculations of Lam and Mahlon \cite{Lam:2000nz},
the nuclear matter density is taken to be uniform.
The nuclear correction factor $\cal{E}$ is determined by
geometry of the nucleus. For cylindrical
and spherical nuclei it is
given by
\begin{equation}
{\cal E}(z)=
\begin{cases}
\displaystyle{\frac{1}{z}} (e^{z}-1) & \text{(cylindrical),} \\ \\
\displaystyle{\frac{3}{z^{3}}} [1-\frac{1}{2}z^{2}+e^{z}(z-1)] &
\text{(spherical).}
\end{cases}
\end{equation}
This function appears in Eq. (\ref{eq:gluon-distro-LM}) as
${\cal E}(v^2 L({\Delta_\perp}))$.
The $v^2$ controls how much the results depend
on the nuclear geometry; it is given by
\begin{equation}
v^2=
\begin{cases}
\displaystyle{ \frac{3 A g^4}{2 \pi R^2} } \approx 24
\pi \alpha_s^2 A^{1/3} R_0^{-2}  & \text{(cylindrical),}
\\ \\
\displaystyle{ \frac{9 A g^4}{4 \pi R^2} } \approx 36
\pi \alpha_s^2 A^{1/3} R_0^{-2} & \text{(spherical).}
\end{cases}
\end{equation}
Here $R=R_0 A^{1/3}$ and we take $R_0=1\,\rm{fm}$.
The difference between the spherical and cylindrical nuclei
is negligibly small. In our calculations, we use the formulae for a cylindrical
nucleus.

Now we are ready to calculate the gluon distribution function
$xg_{A}(x,Q^{2})$
for a nucleus by substituting Eq. (\ref{eq:gluon-distro-LM}) in Eq.
(\ref{eq:gluon-distro}),
\begin{equation}
xg_{A}(x,Q^{2}) = \frac{8 A \alpha_{s}}{\pi^{2}}\int_{\vert
{\boldsymbol{q}_\perp}\vert\leqslant Q}d^{2}
{\boldsymbol{q}_\perp}\int d^{2}
{\boldsymbol{\Delta}_\perp} e^{i {\boldsymbol{q}_\perp}
\cdot {\boldsymbol{\Delta}_\perp}}
{\cal L}(x;{\boldsymbol{\Delta}_\perp}){\cal E}(v^{2}
L({\boldsymbol{\Delta}_\perp})).\label{eq:gluon-distro-final}
\end{equation}
The momentum integration in the equation above is
trivial.
Owing to the longitudinal coordinate $x_{\parallel}$,
the fractional momentum distribution
function $xg_{A}(x,Q^{2})$ in Eq. (\ref{eq:gluon-distro-final})
is x dependent, in contrast
to the results of the original two dimensional MV model.

The region of validity of the 3dMVn model within the
two dimensional parameter space
$(x,Q^{2})$, likewise the original MV model, is restricted
by the assumptions made regarding
the weak coupling limit, color coherence and
color averaging \cite{Lam:2000nz}.

The lower limit of x is set for a gold nucleus by
\begin{equation}
x \sim \frac{A^{-1/3}}{m R_0} \simeq 0.035
\label{lowerlim}
\end{equation}
where $R_0$ is the nucleon size
$\sim 1\,\rm{fm}$. This is the longitudinal momentum
fraction at which the gluons start to resolve the Lorentz
contracted thickness of the
nucleus.

The weak coupling limit imposes an upper
limit on x such that
\begin{equation}
x\lesssim\frac{4}{\pi m R_0}\simeq0.25.
\end{equation}
For very large $Q^2$,
not enough color charge would be seen, hence the
Gaussian approximation for
the color charge fluctuations $\langle\rho^{a}(\boldsymbol{x})
\rho^{b}(\boldsymbol{x}^{\,\prime})\rangle$
would not be valid. Similarly, at the scale of color neutral nucleons
where $Q^{2}$ is very small, again there would not be
enough color charge to average over.
These considerations put limits
on the momentum scale
\begin{equation}
\frac{\pi m }{4 R_0}x\lesssim
Q^{2}\lesssim\frac{4 }{m R_0^3}\frac{1}{x}.
\label{Qvalidity}
\end{equation}
\begin{figure}[t]
\begin{centering}
\includegraphics[scale=0.84]{validity}
\par\end{centering}
\caption{Approximate region of validity of
the 3dMVn model for gold nucleus (A=197). }
\label{figure:validity}
\end{figure}
For $R_0=1\,\rm{fm}$, Eq. (\ref{Qvalidity}) becomes
\begin{equation}
0.16\, x\lesssim Q^{2}\lesssim\frac{0.032}{x}.
\end{equation}

In our calculations, we will take { $Q^{2}=0.55\,{\rm GeV^{2}}$}.
The reason for this choice is that it is at
the upper limit of validity of the 3dMVn model
and at the lower
limit of the parton distribution function to which we will compare.
While the model spans a
larger range in x for $Q^2 < 0.55\,\rm{GeV^2}$
(see Figure \ref{figure:validity}),
no data is available at those scales.
At the energy scale { $Q^{2}=0.55\,{\rm GeV^{2}}$},
x is restricted to the range
\begin{equation}
0.035<x<0.060.
\end{equation}
Figure
\ref{figure:validity} shows the approximate
validity region of the 3dMVn model
for a gold nucleus. We consider gold because
it provides for the
greatest range of validity of the model as
represented by Eq. (\ref{lowerlim}).
The results for lead would be indistinguishable.

The parameters of the 3dMVn model need to be fixed
before performing the numerical integration
in Eq. (\ref{eq:gluon-distro-final}). We take the
nucleon mass $m=0.94\,\mbox{GeV}$,
$R_0=1\,{\rm fm}$ and $A=197$.
The coupling constant $\alpha_{s}$ and the
correlation length $\lambda$
are treated as free parameters.
We evaluate $xg_{A}(x,Q^{2})$ numerically at
$Q^2=0.55\,\rm{GeV^2}$ in the range
$0.035<x<0.060$ for several values of
$\alpha_{s}$ and $\lambda$
in increments of 0.1
in the range
\begin{equation}
 \begin{aligned}
0.1 \leq &\alpha_{s}\leq  1,\\
0.2\,{\rm fm}  \leq & \lambda \leq
2.6\,{\rm fm.}\label{eq:param}
  \end{aligned}
\end{equation}
Then we fit $xg_{A}(x,Q^{2})$ to a form
\begin{equation}
xg_{A}(x,Q^{2})= b \, x^c \, (1-x)^d,
\label{parametrization}
\end{equation}
and find $\{b,c,d\}$ for each set of $\{\alpha_s,\lambda\}$.
After the numerical integration, we obtain
${xg_{197}(x,Q^{2})}$
for the whole nucleus
for several points in the parameter space
$\{\alpha_{S},\lambda \}$.
The scale $Q^{2}=0.55\,{\rm GeV^{2}}$ at
which we calculated $xg_{A}(x,Q^{2})$
is the lower limit of the JR09 data
parametrization for the parton
distribution functions.
For this reason, we compare JR09 with the gluon distribution functions that are calculated and evolved to $Q^{2}=100\,{\rm GeV^{2}}$.

Gluon distribution functions are evolved to different energy scales
with the DGLAP equation
\cite{Gribov:1972ri,Lipatov:1974qm,Dokshitzer:1977sg,Altarelli:1977zs}.
For this purpose, we employ the code
QCDNUM17 \cite{Botje:2010ay}.
We will assume that nucleus is a dilute system of nucleons
and take the nuclear modification factor $R_A=1$
\cite{JalilianMarian:1999cf}. Hence,
\begin{equation}
xg_{p/A}(x,Q^{2})\equiv\frac{1}{A}xg_{A}(x,Q^{2}),
\end{equation}
where the subscript $p/A$ refers to a nucleon in
a nuclear environment.
Henceforth, we will use the gluon distribution function for
a single nucleon in a gold nucleus.
It should be kept in mind that
$g_{p/A}(x,Q^{2})$ is in principle different than a
gluon distribution function for an isolated
proton $g_{p}(x,Q^{2})$ since the nuclear effects
introduce enhancement or shadowing
depending on the energy scale $Q^{2}$.

We run QCDNUM in the variable-flavor number scheme (VFNS)
and at the next-to-next-to-leading order (NNLO($\overline{\text{MS}}$)).
For the DGLAP evolution of the gluons, valence and
sea quark distributions
at the initial scale should also be provided to the evolution
code. At the initial scale $Q^{2}=0.55\,{\rm GeV^{2}}$, the calculated
gluon distribution function is the main input to QCDNUM. As for
the valence and sea quark distributions, the JR09
\cite{JimenezDelgado:2009tv} data
parametrization is utilized. The evolution is
repeated for several
values of $\alpha_{s}$ and $\lambda$
in the range given in Eq. (\ref{eq:param}).

\section{Data Parametrization}

We will compare our results with the parametrization of data by
JR09 \cite{JimenezDelgado:2009tv}.
We use the PDFs that work in the VFNS scheme and at the
NNLO($\overline{\text{MS}}$) order.
The JR09 is valid in the ranges
\begin{equation}
 \begin{aligned}
0.55\lesssim & Q^{2}\lesssim10^{8}\,{\rm GeV^{2}}, \\
10^{-9}\lesssim & x\lesssim1.
  \end{aligned}
\end{equation}

The JR09 provides
PDFs for separate nucleons $xg_{p}(x,Q^{2})$ but
it does not contain any information about the nuclear modification
function $R_A$. Our results, on the other hand, are for a nucleon in
a nuclear environment $xg_{p/A}(x,Q^{2})$; therefore,
$R_A$ should be, in principle, taken into account. However,
the nuclear PDFs (for example nCTEQ
\cite{Schienbein:2009kk}) are available only for
$Q^{2}>1\,{\rm GeV^{2}}$.
As pointed out earlier, since the DGLAP equation mixes quarks and
gluons during the evolution, we need the valence and sea quark distribution
functions at the initial energy scale so that the gluon distribution function
calculated from the model can be evolved to the higher energies.
Hence, for consistency, we will utilize
the quark PDFs provided by JR09 at $Q^{2}=0.55\,{\rm GeV^{2}}$ as an
input to the DGLAP evolution.
The discrepancy
between JR09 and nCTEQ \cite{Schienbein:2009kk},
encoded in the nuclear modification factor,
is at most $5\%$ at $Q^{2}=25\,{\rm GeV^{2}}$.
At $Q^{2}=100$ and $1000\,{\rm GeV^{2}}$ the discrepancy is
negligibly small. For that reason, we think
the lack of information regarding
$R_A$ will not cause a significant error in our analysis.

\section{Best Fit Parameters}
Now we seek the sets of parameters $\{\alpha_s,\lambda \}$
that produce the best fits. For this reason,
we compare the gluon distribution functions,
calculated and evolved to $Q^{2}=100\,{\rm GeV^{2}}$,
with JR09 at the same energy.

In the DGLAP evolution to higher energies, the contribution
to the radiation at a particular value of x always
comes from the sources from the
larger x region.
Hence, the effect of gluons calculated from 3dMVn
in the range $0.035<x<0.060$ at
$Q^{2}=0.55\,{\rm GeV^{2}}$ will be more prominent
at smaller x values at
$Q^{2}=100\,{\rm GeV^{2}}$. For this purpose,
we will compare the evolved model with JR09 in the range
\begin{equation}
0.015<x<0.04.
\label{fitregion}
\end{equation}
The gluons from $x<0.015$ do not enter the DGLAP evolution
and hence they do not contribute to the region in Eq. (\ref{fitregion}).
The gluons from $0.015<x<1$ will enter the
DGLAP evolution. We argue that the
parametrization in Eq. (\ref{parametrization}) is a
good approximation
for the large-x gluons.

Our analysis shows that the best fit occurs at the values
$\alpha_{s}=0.5$ and
$\lambda=1.8\,{\rm fm}$. In Figure \ref{fig:modelanddata}
we show a comparison between the JR09 data parametrization
and 3dMVn model at $Q^{2}=100\,{\rm GeV^{2}}$ for
$\alpha_{s}=0.5$ and
$\lambda=1.8\,{\rm fm}$.
We find almost identical plots at other $Q^2$
values in the range
${0.015<x<0.04}$.
\begin{figure}[t]
\begin{centering}
\includegraphics[scale=0.7]{modelanddata.eps}
\par\end{centering}
\caption{Comparison of the gluon distribution
function of a nucleon from the 3dMVn model
for various values of $\{\alpha_s,\lambda\}$
with the JR09 data parametrization at
$Q^2= 100\,\rm{GeV^2}$. The horizontal axis is
in logarithmic scale.
For distances larger than the correlation
length $\lambda$, gluons are not
correlated and hence the color neutrality condition
is satisfied.
The model is reliable in the
fit region bounded by the two vertical lines.
The uncertainty in JR09 is shown with
an error band and it is about $5\%$.
The discrepancy between the best fit curve
($\alpha_s=0.5$ and $\lambda=1.8\,\rm{fm}$)
and the JR09 in the fit region is only $2\%$.
}
\label{fig:modelanddata}
\end{figure}

At other values of $\alpha_s$ and $\lambda$,
we find that the model underestimates the data.
As $\alpha_s$ increases, the 3dMVn curve increases
and the disagreement between the model and JR09 decreases.
For values $\alpha_s \geq 0.5$
the model does not change significantly.
Hence, we take
${\alpha_s=0.5}$ since this is the smallest value of $\alpha_s$
for which a good fit can be obtained.
Interestingly, this ``freezing'' behavior of the
strong coupling constant
particularly at $\alpha_s=0.5$,
imposed in other approaches by hand
(see \cite{Dumitru:2011wq,Albacete:2012xq}),
arises naturally in the 3dMVn model.
Adjusting $\lambda$ changes the model curve only slightly.
For $\alpha_s=0.4$ and
values of $\lambda$ greater than $1.8\,\rm{fm}$,
we do not find
a fit as good as the one for $\alpha_s=0.5$
and $\lambda=1.8\,\rm{fm}$.

There are various sources of uncertainty
in the MV model which are shared with the 3dMVn model.
In the paradigm of strong classical color fields,
one wishes to employ solutions of
the classical Yang-Mills equation which is given
in covariant gauge by
\begin{equation}
(\boldsymbol{\nabla}_{\perp}^2+ \partial_{\parallel}^2)A^\nu=
g J^\nu
+
2ig[A^\mu,\partial_\mu A^\nu]
-ig[A_\mu,\partial^\nu A^\mu]
+
g^2[A_\mu,[A^\mu,A^\nu]]
,
\label{cym}
\end{equation}
where $\partial_{\parallel}^2=-2 \partial_{+} \partial_{-}$.
The second and third terms
on the right hand side of the Eq. (\ref{cym}) are
responsible for the
three gluon interaction vertex, and the last term
is for the four gluon interaction vertex in the quantized
theory. In the framework of Color Glass Condensate,
only the linearized version of Eq. (\ref{cym}) is used.
The nonabelian feature of the model sets in when
the solutions of the linearized equation in the covariant
gauge are transformed to the light-cone gauge by
means of the full gauge transformation including the
nonabelian term \cite{Lam:1999wu}
\begin{equation}
A_{\rm{LC}}^{\mu}(x)={\rm{U}}(x)
A_{\rm{cov}}^{\mu} {\rm{U}}^{-1}(x)
- \frac{i}{g} [\partial^{\mu}  {\rm{U}}(x)] {\rm{U}}^{-1}(x).
\end{equation}
The reason why the nonlinear terms in Eq. (\ref{cym})
are omitted
is simply because the solutions of the fully nonlinear theory
are not known. In addition, the Green's function for gluons,
which is essential in the calculation of
{$\langle A_{i}^{a}(\boldsymbol{x})A_{i}^{b}(
\boldsymbol{x}^{\,\prime})\rangle$} from
{$\langle\rho^{a}(\boldsymbol{x})\rho^{b}(
\boldsymbol{x}^{\,\prime})\rangle$},
can only be defined in the linear theory.
The gluon propagator is calculated from
\begin{equation}
(\boldsymbol{\nabla}_{\perp}^2+ \partial_{\parallel}^2)
G(\boldsymbol{x})= \delta^3 (\boldsymbol{x}).
\end{equation}

The origin of the infrared divergences in the MV model
lies here and this is also why confinement effects need to be
introduced by hand by using colored noise. In principle,
the correlation function with colored noise
in Eq. (\ref{eq:Yukawa})
should be calculable from the fully nonlinear theory.

The linear
approximation can only be justified if $g \ll 1$
so that
the nonlinear terms can be neglected.
However,
the strong coupling constant is not expected to
be so small and nonlinear terms are as
important as the source term.
Hence, the linear approximation
of Eq. (\ref{cym}) is the main source of the uncertainty
in any realistic analysis.

Another source of uncertainty
may be the assumption that
the nucleus is taken to be much larger
than a nucleon, $R \gg R_0$. For a gold
nucleus, the corrections may be as large as $17\%$ since
$R_0/R \sim A^{-1/3} \sim 0.17. $
Also quantum corrections at
smaller values of x may be important
even in the weak coupling limit. These corrections
are of the form $\alpha_s \ln (1/x)$ and
discussed in the references
 \cite{Ayala:1995kg,Ayala:1995hx,Lam:2000nz}.
Lastly, although we estimate that the error due to
neglecting the nuclear effects should be small,
it would be worth of repeating the analysis
presented in this paper by using nuclear PDFs
once they are available for $Q^2 \sim 0.55\, \rm{GeV^2}$.

\section{Summary and Outlook}

The framework of Color Glass Condensate allows
one to calculate the vector field correlation
function. It can be used to calculate the gluon
distribution function of the ultrarelativistic nuclei and
the initial energy density distribution due to the
interacting classical color fields
produced by the colliding nuclei in heavy ion collisions.

We have examined a three dimensional color neutral version of
the McLerran-Venugopalan model. The 3dMVn
model is finite in the infrared
and therefore an infrared cutoff is not needed.
In addition, the results of this model are
x dependent due to the intrinsic three dimensional
treatment of the nucleus
in contrast to the approximation of infinitely thin nucleus.

In order to explore the parameter space of these
two variables, we have calculated the gluon distribution
function for several values of $\alpha_s$ and $\lambda$.
The originality of this work is to compare our calculations directly with parametrization of the data to determine the free parameters of the model. We have found the best fit between
the gluon distribution function from the JR09 parametrization
and the one calculated from the 3dMVn model occurs at
$\alpha_s=0.5$ and $\lambda=1.8\,\rm{fm}.$
At other values of these two parameters the model underestimates
the data. This may be due to the uncertainty
in the assumptions and
rough estimates made during the construction of the model
as well as the other uncertainties discussed in the previous
section.
We have also found that the 3dMVn model had an intrinsic freezing behavior that the gluon distribution function froze at $\alpha_s=0.5$ and remained unchanged for $\alpha_s>0.5$.

The assumption that the color charge is normally distributed throughout the nucleus lies at the heart of the MV and 3dMVn models. Besides the normal (Gaussian) distribution of the color charge, the nucleonic inner structure of the nucleus is implemented via colored noise by introducing correlations for distances smaller than a nucleon size. In this manner, the fluctuations in the positions of nucleons in a nucleus have not been treated separately from the dynamical color charge fluctuations. In other words, the effect of confinement is realized through short range correlations in $\langle\rho^{a}(\boldsymbol{x})\rho^{b}(\boldsymbol{x}^{\,\prime})\rangle$ by colored noise rather than considering the nucleus as a collection of individual nucleons. An alternative to this might be an event-by-event Monte Carlo sampling of the distribution of nucleons \cite{Schenke:2010rr,Schenke:2011tv,Schenke:2012wb,Schenke:2012hg}. For a finite nucleus, the results from event-by-event fluctuations in positions of the sampled nucleons may differ from the analytical calculations presented in this paper. For a very large nucleus, we expect the results from both approaches to agree.

The final product of this work is the unintegrated gluon distribution (also known as the ``gluon number density'') given in Eq. (\ref{eq:gluon-distro-LM}) of which parameters are fixed in the previous section. This three dimensional quantity is the main ingredient for various applications of the CGC concept. It can be picked up by the practitioners of CGC and be used immediately in realistic calculations other than calculation of the gluon distribution function.

In this work, we calculated the gluon distribution function since it was a quantity that could be easily compared with data.
However,
the ultimate goal was not
calculating the gluon distribution function per se,
but to attain a working model for further use,
particularly in the calculation of the initial energy density
distribution in heavy ion collisions.
This work is in progress.
In the future work, we plan to calculate the initial
energy density in heavy ion collision in the framework of the 3dMVn model.

\section*{Acknowledgments}
We thank Joseph Kapusta and Rainer J. Fries for
useful discussions and their critical comments on the manuscript.
We also thank Michiel Botje for correspondence regarding
his QCDNUM code. This work was supported by the U.S. DOE Award
No. DE-FG02-87ER40328.

\raggedright


\section{\label{sec:intro} Introduction}

  The idea of breaking the inversion symmetry by some complex magnetic
order has attracted a great deal of attention. It gives rise to the phenomenon of
improper multiferroicity, where the ferroelectric (FE) polarization is
induced by the magnetic order and, therefore, can be controlled by the magnetic field.
Alternatively, one can change the magnetic structure by applying the electric field.
This makes multiferroic materials promising for the creation of the new generation of
electronic devices.\cite{MF_review,spiral_review}

  One of magnetic textures, that breaks the inversion symmetry, is the spin spiral.
Moreover, this property is universal, in a sense that it can take place in all types
of magnetic compounds, irrespectively of their symmetry, as long as the
the spin-spiral order is established.\cite{Sandratskii} Probably, this is the reason why
the idea of the spin-spiral alignment became so popular in the field of multiferroics:
many properties of such compounds are interpreted from the viewpoint of
spin-spiral order,\cite{KNB,Mostovoy,Sergienko} and the search for new multiferroic
materials is frequently conducted around those with the spin-spiral texture.\cite{spiral_review}

  Nevertheless, the spin-spiral alignment alone is insufficient for having a FE activity:
although the spin-spiral order breaks the inversion symmetry, one can always find
some appropriate uniform rotation of spins, which, after combining with the spatial inversion,
will transform the spin-spiral texture to itself (see Appendix~\ref{sec:appendix}). The
FE polarization in such a situation will be equal to zero. In order to make it finite
(and in order to fully break the inversion symmetry), one should ban the possibility of uniform rotations
in the system of spins.
This can be done by the relativistic spin-orbit (SO) interaction, and today it is commonly
accepted in the theory of multiferroic materials that the spin-spiral alignment
should be always supplemented with the SO coupling.\cite{KNB,Sergienko}
However, as long as the SO coupling is involved, the spin-spiral texture becomes deformed.
Then, we face the question whether the FE activity is caused by the spin-spiral order itself
or by the deviation from it.\cite{antispiral}
The answer to this question is very important for understanding the origin of multiferroicity.

  MnWO$_4$ is an important material in the field of multiferroics.
Although the transition temperature
to the multiferroic phase is low and
the FE polarization is weak, the material is fundamentally important, because: (i) its multiferroic
phase has a spin-spiral texture (or is believed to be the spin spiral);
(ii) the relative directions of the spin-spiral propagation vector,
the spin rotation axis, and the FE polarization seem to be consistent with predictions based
on the spin-spiral theory.\cite{Taniguchi,Arkenbout,Heyer} Thus, MnWO$_4$ is regarded as one of
successful manifestations of the spin-spiral theory and, actually, the discovery of the FE activity in MnWO$_4$
was guided by this theory.

  More specifically, MnWO$_4$ crystallizes in the monoclinic $P2/c$ structure.
It exhibits three successive antiferromagnetic (AFM)
transitions: at $T_{\rm N3}$$\sim 13.5$ K, $T_{\rm N2}$$\sim 12.5$ K, and $T_{\rm N1}$$\sim 7.6$ K.
The AF1 phase, realized below $T_{\rm N1}$, is nearly collinear and the directions of spins
alternate along the monoclinic ${\bf b}$ axis as $\uparrow \uparrow \downarrow \downarrow$
(the reason why the AF1 phase is also called as the `$\uparrow \uparrow \downarrow \downarrow$ phase').
The AF2 phase ($T_{\rm N1}$$<T$$<$$T_{\rm N2}$) is incommensurate
and is typically ascribed to the spin spiral with the
propagation vector
${\bf q}_{\rm AF2} = (-$$0.214,1/2,0.457)$ (measured in units of reciprocal lattice translations).
The spin moments lie in the plane formed by the easy magnetization direction in the plane
$\boldsymbol{ac}$ and the $\boldsymbol{b}$ axis.
The AF3 phase ($T_{\rm N2}$$<$$T$$<$$T_{\rm N3}$) is characterized by a collinear sinusoidally
modulated AFM order with the same ${\bf q}$.
The AF2 phase is ferroelectric. The polarization vector is
parallel to the monoclinic $\boldsymbol{b}$ axis: ${\bf P} = (0,P_b,0)$.
$P_b$ is relatively small: it takes a maximal value ($\sim$ 50 $\mu$C/m$^2$) at around $T_{\rm N1}$, and
then monotonously decreases and vanishes at around $T_{\rm N2}$.\cite{Taniguchi}
Other phases display no sign of FE activity.

  In this work, we will argue that the multiferroicity in MnWO$_4$ has a more complex origin. Namely, we will show that:
(i) It is not sufficient to have a simple spin-spiral texture; (ii) The spin spiral should be deformed by some
conflicting interactions, existing in the system. In MnWO$_4$, these are Dzyaloshinskii-Moriya (DM) interactions
and isotropic exchange interactions. The conflict of these two interactions breaks the inversion
symmetry and gives rise to the FE activity in the AF2 phase.

\section{\label{sec:method} Method}

  We follow the same strategy as in our previous publications,
devoted to multiferroic manganites.\cite{antispiral} Our basic idea is to construct a realistic
low-energy (Hubbard-type) model for the Mn $3d$ bands, which would be sufficient for describing basic
magnetic properties of MnWO$_4$ at a semi-quantitative level. The model is constructed in the basis of
Wannier orbitals, starting from the electronic structure calculations
in the local density approximation (LDA) with and without the relativistic SO interaction. The details
can be found in the review article (Ref.~\onlinecite{review2008}).
All calculations have been performed using experimental parameters of the
crystal structure.\cite{Lautenschlaeger}
Since heavy W atoms is an important source of relativistic effects,
the SO coupling was included already at the level of ordinary LDA calculations, and thus
obtained electronic structure was used as the starting point for the
construction of the low-energy model.
Thus, although the W sites were not explicitly included to the model, the relativistic effects,
associated with these sites, were allowed to contribute to the dispersion
of the Mn $3d$ bands and, in this way, were taken into consideration during the analysis of
the low-energy model. All obtained parameters of the model Hamiltonian
are collected in Ref.~\onlinecite{SM}. After the construction, the model was solved in the
mean-field Hartee-Fock (HF) approximation.

  We would also like to briefly comment on merits and demerits of our technique in comparison with
the first-principles approach:

  (i) It is certainly true that the construction of the model Hamiltonian is based on some
additional approximations,\cite{review2008} and is sometimes regarded as a step back
in comparison with first-principles calculations. On the other hand, it would not be right to think
that the first-principles calculations for the transition-metal oxides are free of any approximations.
The necessity to treat the problem of on-site Coulomb correlations, which is frequently formulated
in terms of the LDA$+$$U$ approach (with some phenomenological correction ``$+$$U$'', borrowed from the
Hubbard model) make it similar to the model approach: it also relies on additional approximations,
although of a different type. The typical approximations are the
choice of the parameter $U$ and the form of the double-counting term, which are both
ill-defined in LDA$+$$U$.\cite{PRB98} Thus, at the present stage, it would not be right
to say that one techniques is more superior than other: rather, they provide a complementary information
for the analysis of material properties of the transition-metal oxides.
As we will see below, results of our model analysis are well consistent with available
first-principles calculations for MnWO$_4$.\cite{Tian,Shanavas}

  (ii) The search of the ground state for noncollinear magnetic textures can be very
time consuming, even at the level of mean-field HF approximation: the rotations of magnetic moments
towards new equilibrium positions after including the SO interactions can be very slow and require tens of thousands of iterations,
as in the case of MnWO$_4$. In such a situation, the model Hamiltonian approach is very useful,
because it allows us to find a fully self-consistent
magnetic structure, which is not always accessible for
the ordinary first-principles calculations.\cite{antispiral}

 (iii) The model Hamiltonian approach can be very useful for the analysis and interpretation of results
of complex electronic structure calculations and the experimental data. Particularly, in this work it will
help us to elucidate the microscopic origin of the magnetic inversion symmetry breaking in MnWO$_4$.
We will also show that the experimental behavior of long-range magnetic interactions
in MnWO$_4$ (Ref.~\onlinecite{Ye}) can be naturally understood in the framework of the superexchange (SE)
theory and reflects similar behavior of the transfer integrals, derived for the low-energy model.

\section{\label{sec:results} Results and discussions}

  First, we solve the low-energy model in the HF approximation,
by assuming the collinear ferromagnetic (FM) alignment of spins, and drive parameters of
isotropic exchange interactions. For these purposes we employ the theory of
infinitesimal spin rotations.\cite{review2008,Liechtenstein}
The procedure corresponds to the local mapping of the change of the one-electron energy onto the spin Hamiltonian
of the form
${\cal H}_S = -\sum_{\langle ij \rangle} J_{ij} {\bf e}_i \cdot {\bf e}_j$, where
${\bf e}_i$ is the \textit{direction} of spin at the site $i$ and the
summation run over \textit{inequivalent pairs} of sites.
The results of these calculations are explained in Fig.~\ref{fig.J}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=7.1cm]{figure1a.eps}
\includegraphics[width=5.8cm]{figure1b.eps}
\end{center}
\caption{\label{fig.J}
(Color online) (a) Lattice of Mn sites with the notations of
isotropic exchange interactions
(numbered in the increasing order of interatomic distances).
Note that the $P2/c$ structure of
MnWO$_4$ has two Mn sublattices, which are shown by different colors
and denoted as `$I$' and `$II$', respectively. These sublattices are
transformed to each other by the inversion operation.
(b) Distance-dependence of isotropic exchange interactions:
results of calculations, using the theory of infinitesimal
spin rotations near the ferromagnetic state (denoted as `inf') and the
theory of superexchange interactions with $\Delta_{\rm ex} = 5$ eV (denoted as `SE'), in comparison with
the experimental data from Ref.~\protect\onlinecite{Ye}.
}
\end{figure}
Alternatively, one can apply the theory of superexchange (SE) interactions,
by considering the energy gain caused by virtual hoppings ($t_{ij}^{mm'}$) in
the second order of perturbation theory.\cite{PWA}
For the $d^5$ configuration of Mn$^{2+}$,
this expression is extremely simple:
$J_{ij} = -$$\sum_{mm'} t_{ij}^{mm'} t_{ji}^{m'm}/\Delta_{\rm ex}$, where $\Delta_{\rm ex}$
is the intraatomic splitting between the majority- and minority-spin states and the summation runs
over all $3d$ orbitals: $m$ $(m') =$ $xy$, $yz$, $3z^2$$-$$r^2$, $zx$, and $x^2$$-$$y^2$.
Thus, all SE interactions are expected to be antiferromagnetic.
Then, by taking $\Delta_{\rm ex} = 5$ eV and using the values of
$t_{ij}^{mm'}$ obtained for the low-energy model (Ref.~\onlinecite{SM}), one can find that the SE theory
excellently reproduces the
results of the more general theory of infinitesimal spin rotations.
Moreover, the value of $\Delta_{\rm ex} = 5$ eV
is well consistent with typical estimates in the framework of the Hubbard model, where
$\Delta_{\rm ex} \approx U$$+$$4J_{\rm H}$ and the parameters of averaged on-site Coulomb repulsion
$U$ and the exchange interaction $J_{\rm H}$,
derived for the low-energy model, are $1.8$ and $0.8$ eV, respectively.\cite{SM}
Thus, all estimates are consistent with each other. Therefore, we may conclude that:
(i) the physically relevant mechanism, responsible for isotropic exchange interactions
in MnWO$_4$, is the superexchange; and (ii) the distance-dependence of
$J_{ij}$ reflects the behavior of the transfer integrals $t_{ij}^{mm'}$.

  Our calculations capture main details of the experimental
distance-dependence of $J_{ij}$ (Fig.~\ref{fig.J}).\cite{Ye} This is a very important finding, which means that
in our model we should be able to reproduce the correct
magnetic structures of MnWO$_4$.
The negative aspect is that our parameters $J_{ij}$ are `too antiferromagnetic' (AFM).
For example, using the values of $J_{ij}$ reported in Fig.~\ref{fig.J},
the Curie-Weiss temperature can be estimated as
$\theta_{\rm CW} \approx \sum_j J_{ij}/3 k_{\rm B}$.
It yields $\theta_{\rm CW}$$= -$$265$ K, which exceeds the experimental value
$-$$78$ K (Ref.~\onlinecite{Taniguchi}) by factor three.
The discrepancy cannot be resolved simply by changing the value of the on-site Coulomb repulsion $U$,
which is frequently regarded as an ill-defined parameter. Particularly, in order to explain the existence of
FM interactions, which were observed experimentally in some of the bonds,\cite{Ye} we need an additional
mechanism on the top of our model. Such a mechanism can be related to the
magnetic polarization of the oxygen band, similar to orthorhombic manganites.\cite{JPSJ}
The first-principles GGA$+$$U$ calculations
(where GGA stands for the generalized gradient approximation -- an extension of LDA)
also substantially overestimate $| J_{ij} |$ and $| \theta_{\rm CW} |$.\cite{Tian}
Although these calculations take into account the effect of the oxygen band,
the disagreement is probably caused by additional approximations for the double-counting energy,
which is ill-defined in GGA$+$$U$.\cite{PRB98} Finally, the experimental parameters $J_{ij}$
themselves may be sensitive to the details of fitting of the spin-wave dispersion. In this respect,
we would like to mention that earlier experimental data were interpreted
in terms of rather different set of parameters $J_{ij}$.\cite{Ehrenberg}
Thus, at the present stage, taking into account the complexity of
the problem of interatomic magnetic interactions in MnWO$_4$,
the agreement with the experimental data in Fig.~\ref{fig.J} can be regarded as satisfactory.

  Then, we search for the magnetic ground state of MnWO$_4$ without relativistic SO coupling.
For these purposes we employ the spin-spiral formalism, which is based on the generalized Bloch
theorem.\cite{Sandratskii} The idea is to combine the lattice translations with
rotations in the spin subspace. It allows us to treat the incommensurate magnetic textures
with an arbitrary vector ${\bf q}=(q_a,1/2,q_c)$ in the reciprocal space
such as if the magnetic primitive cell was of the same size as the
crystallographic one.\cite{Sandratskii} The results of these
calculations, which were also performed in the HF approximation
for the low-energy model, are summarized in Fig.~\ref{fig.Eq}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=10cm]{figure2.eps}
\end{center}
\caption{\label{fig.Eq}
(Color online) Results of Hartree-Fock calculations
for the spin-spiral states with ${\bf q} = (q_a,1/2,q_c)$: total energies versus $q_a$ for different values of $q_c$
(in units of $2\pi/c$).
The energy of the collinear $\uparrow \uparrow \downarrow \downarrow$ state is marked by the symbol $\otimes$.
}
\end{figure}
For ${\bf q} = (-$$1/4,1/2,1/4)$, the
collinear $\uparrow \uparrow \downarrow \downarrow$ and noncollinear spin-spiral states
are nearly degenerate. Nevertheless, the total energy of the spin spiral continues to decrease
with the increase of $q_a$. Therefore, the `nonrelativistic' ground state of MnWO$_4$ is expected to be
the spin spiral with ${\bf q}_{\rm NR} \approx (-$$0.175,1/2,0.45)$, which is lower in energy than the
$\uparrow \uparrow \downarrow \downarrow$ state
by about $0.5$ meV/Mn.
Moreover, none of these states is ferroelectric: although the spin spiral breaks
the inversion operation, the latter can be always combined with an appropriate spin rotation, which
will transform
the inverted spin spiral to the original one. Thus,
the FE polarization will vanish (see Appendix~\ref{sec:appendix}).

  Therefore, we should find the answer to the following questions: (i) Which interaction stabilizes
the AF1 structure and make it the ground state of MnWO$_4$
at low $T$? (ii) What is the mechanism of the inversion-symmetry breaking in
the AF2 state, which yields finite $P_b$?
For these purposes we consider the relativistic SO interaction.
It gives rise to such important ingredients as the single-ion anisotropy and DM interactions.

  The effect of single-ion anisotropy on the magnetic texture of MnWO$_4$ can be studied by
enforcing the atomic limit and setting all transfer integrals equal to zero.\cite{PRB12}
Generally, the single-ion anisotropy is expected to be small for the nearly spherical
$d^5$ configuration of Mn$^{2+}$.
Since the symmetry operation $\{ C_b^2|\boldsymbol{c}/2 \}$ of the $P2/c$ space group transforms each Mn site to itself
(here, $C_b^2$ stands for the $180^\circ$ rotation around the $\boldsymbol{b}$ axis and $\boldsymbol{c}/2$
is the additional translation),
the single-ion anisotropy will tend to align the spins either parallel
to the $\boldsymbol{b}$ axis ($|| \boldsymbol{b}$)
or in the $\boldsymbol{ac}$ plane ($\in \boldsymbol{ac}$). In the latter case, the spins are canted
off the $\boldsymbol{a}$ axis by about $41^\circ$,
which is close to the experimental value of $35^\circ$.\cite{Taniguchi}
The easy magnetization direction corresponds to the in-plane configuration. However, the energy difference between
hard ($|| \boldsymbol{b}$) and easy ($\in \boldsymbol{ac}$) magnetization
directions is small (less that $0.05$ meV/Mn), so that in the spin-spiral texture it
can be overcome by the energy gain caused by the isotropic exchange interactions.

  The $P2/c$ structure of MnWO$_4$ has two Mn-sublattices (denoted as `$I$' and `$II$' in Fig.~\ref{fig.J}),
which are transformed to each other by the inversion operation.
Therefore, the DM interactions will obey the following symmetry rules:\cite{DM}
\begin{itemize}
\item[(i)]
all interactions between different sublattices are equal to zero;
\item[(ii)]
the DM interactions in the sublattices `$I$' and `$II$' are related by the identity
${\bf d}^{I}_{\boldsymbol{R}} = {\bf d}^{II}_{-\boldsymbol{R}} = -$${\bf d}^{II}_{\boldsymbol{R}}$,
for any translation $\boldsymbol{R}$, connecting the sites within one Mn sublattice
(note that ${\bf d}$ is the axial vector).
\end{itemize}
From now on, it is more convenient to use the
extended notations and recall that each site of the lattice $i$ is specified by its
position $\boldsymbol{\tau}$ in the primitive cell and the translation $\boldsymbol{R}$.

  Moreover, due to the symmetry operation $\{ C_b^2|\boldsymbol{c}/2 \}$,
all DM vectors lie in the $\boldsymbol{ac}$ plane, i.e., similar to the easy magnetization
direction, obtained from the analysis of
the single-ion anisotropy.
The parameters of DM interactions can be evaluated by considering the mixed type of the perturbation theory
expansion with respect to the SO coupling and infinitesimal rotations of
spins near the FM state with subsequent mapping of obtained results for the change of the
one-electron energy onto the spin Hamiltonian of the form
$\sum_{\langle \boldsymbol{R} \boldsymbol{R}' \rangle } {\bf d}_{\boldsymbol{R} - \boldsymbol{R}'}
[{\bf e}_{\boldsymbol{R}} \times {\bf e}_{\boldsymbol{R}'} ]$, constructed
for each magnetic sublattice.\cite{PRL96} It yields the following parameters of DM interactions (in meV):
${\bf d}^{I}_{\boldsymbol{a}} = (-$$0.01,0,0.01)$, ${\bf d}^{I}_{\boldsymbol{c}} = (-$$0.02,0,0.02)$,
and ${\bf d}^{I}_{\boldsymbol{a}+\boldsymbol{c}} = (0.01,0,-$$0.01)$. Interactions in other bonds are
considerably smaller.

  The equilibrium magnetic structure, obtained
after switching on the relativistic
SO interaction
in the collinear
AF1 state, is explained in Fig.~\ref{fig.uudd}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=10cm]{figure3.eps}
\end{center}
\caption{\label{fig.uudd}
(Color online) (a) The
magnetic structure AF1, as obtained in the mean-field HF calculations for the
low-energy model with relativistic SO interaction, and (b) its projection onto the $\boldsymbol{ab}$ plane.
}
\end{figure}
Due to the single-ion anisotropy,
all spins are
confined mainly in the $\boldsymbol{ac}$ plane and
canted off the $\boldsymbol{a}$ axis by about $46^\circ$
(which is close to $41^\circ$, obtained in the atomic limit).
The DM interactions yield an additional canting
out of the $\boldsymbol{ac}$ plane. The corresponding rotational force ${\bf f}$, acting at
some Mn spin from its neighboring sites in the same
sublattice, is given by
${\bf f} = \sum_{\boldsymbol{R}} [ {\bf d}_{\boldsymbol{R}} \times {\bf e}_{\boldsymbol{R}} ]$.
Since ${\bf d}_{-\boldsymbol{R}} = -$${\bf d}_{\boldsymbol{R}}$,
in order to contribute to ${\bf f}$, the directions of neighboring spins
should also satisfy the condition ${\bf e}_{-\boldsymbol{R}} = -$${\bf e}_{\boldsymbol{R}}$.
In the $\uparrow \uparrow \downarrow \downarrow$ texture, such a situation takes place for
$\boldsymbol{R} =$ $\boldsymbol{a}$ and
$\boldsymbol{a}$$+$$\boldsymbol{c}$, but not for
$\boldsymbol{R} =$ $\boldsymbol{c}$.
Moreover, since both ${\bf d}_{\boldsymbol{R}}$ and
${\bf e}_{\boldsymbol{R}}$ lie in the $\boldsymbol{ac}$ plane, the force
${\bf f}$ is parallel to the $\boldsymbol{b}$ axis.
Finally, due to the alternation of ${\bf e}_{\boldsymbol{R}}$
in the $\uparrow \uparrow \downarrow \downarrow$ texture,
the directions of such forces will also alternate. All in all, this explains fine details of the
magnetic structure in Fig.~\ref{fig.uudd}b.
Very importantly, the isotropic exchange interactions, single-ion anisotropy, and DM interactions
in the AF1 phase
do not conflict with each other in a sense that all of them,
even jointly, continue to respect the crystallographic $P2/c$ symmetry.
Therefore, the inversion symmetry is preserved and the polarization is equal to zero, even
after including the SO interaction.

  The effect of SO interaction on the spins-spiral alignment can be best understood for
${\bf q} = (-$$1/4,1/2,1/2)$ (see Fig.~\ref{fig.spiral4}).
\begin{figure}[h!]
\begin{center}
\includegraphics[width=15cm]{figure4.eps}
\end{center}
\caption{\label{fig.spiral4}
(Color online) (a) The deformed spin-spiral
texture, as obtained in the mean-field HF calculations for the
low-energy model with relativistic SO interaction for
${\bf q} = (-$$1/4,1/2,1/2)$, and (b) its projection onto the $\boldsymbol{ac}$ plane.
}
\end{figure}
In this structure,
the easy magnetization direction (${\bf e} \in \boldsymbol{ac}$) at certain
Mn-site alternates with the hard magnetization direction (${\bf e} || \boldsymbol{b}$)
at its neighboring sites. Thus, there is a
conflict between isotropic exchange interactions and the single-ion anisotropy. However, since the latter is small,
this magnetic structure can be realized at a low energy cost.
Then, due to the DM interactions,
each site with ${\bf e} || \boldsymbol{b}$ will experience a force from its
neighboring sites with ${\bf e} \in \boldsymbol{ac}$. Since ${\bf d} \in \boldsymbol{ac}$,
this force will be parallel to the
$\boldsymbol{b}$ axis. Alternatively, each site with ${\bf e} \in \boldsymbol{ac}$
will experience the rotational force from its neighboring sites with ${\bf e} || \boldsymbol{b}$.
Since ${\bf d} \in \boldsymbol{ac}$,
this force will also lie in the $\boldsymbol{ac}$ plane.

  Furthermore, since
${\bf d}^{I}_{\boldsymbol{R}} = -$${\bf d}^{II}_{\boldsymbol{R}}$,
the forces in different sublattices will act
\textit{in the opposite directions}.
This means that, if some DM vector ${\bf d}^{I}_{\boldsymbol{R}}$ in the first sublattice
will tend to support the right-handed spiral, its counterpart in the second sublattice
(${\bf d}^{II}_{\boldsymbol{R}}$) will support the left-handed spiral.
However, the spin-spiral texture, that is formed by the isotropic exchange interactions
$J_{ij}$, is always either right- or left-handed \textit{in both magnetic sublattices}.
Therefore, we have a conflict between isotropic and DM interactions: if in one magnetic
sublattice, the effect of these two interactions is added, in other magnetic
sublattice, it will be subtracted. All these tendencies are clearly seen
in fine details of the magnetic structure, obtained in the HF calculations with SO coupling (Fig.~\ref{fig.spiral4}).
First, the canting of spins, lying in the $\boldsymbol{ac}$ plane, appears to
be different for
the sublattices $I$ and $II$:
$21^\circ$ and $34^\circ$, respectively, relative to the $\boldsymbol{a}$ axis.
Second, as discussed above, the spins with ${\bf e} || \boldsymbol{b}$ will experience
the force ${\bf f}$, which is also parallel to $\boldsymbol{b}$. However,
since ${\bf d}^{I}_{\boldsymbol{R}} = -$${\bf d}^{II}_{\boldsymbol{R}}$,
the vectors ${\bf e}$ and ${\bf f}$ at the same site will be either parallel or antiparallel,
depending on the magnetic sublattice. Thus, the magnitude
of spin
magnetic moment with ${\bf e} || \boldsymbol{b}$
will also depend on the magnetic sublattice. Of course, this change of the
magnetic moments is small (about $0.01$ \%, according to our HF calculations).
Nevertheless, it does take place and is just another manifestation
of the conflict between isotropic and DM
interactions.

  Thus, the conflict between isotropic and DM interactions in the AF2 phase
makes the sublattices $I$ and $II$ inequivalent. Since in the
$P2/c$ structure, these two sublattices were connected
by the inversion operation,
their inequivalency means that the inversion symmetry is broken and the actual symmetry of the
AF2 phase is lower than $P2/c$.
Furthermore, the homogeneous spin-spiral alignment itself is
deformed by the DM interactions and,
strictly speaking, the AF2 phase is \textit{no longer the spin spiral}. The magnitude of this
deformation can be seen in the inset of Fig.~\ref{fig.spiral4}.

  Then, we discuss relative stability of different magnetic phases
with the SO interaction. In this case, the generalized Bloch
theorem is not applicable and there only possibility is to work with the
supercell geometry.\cite{Sandratskii} Therefore, we were able to consider
only solutions with ${\bf q} = (q_a,1/2,1/2)$ and the rational $q_a$$=$
$-$$1/3$, $-$$1/4$, and $-$$1/5$.
The results can be summarized as follows. Amongst AF2 states, the one
with $q_a$$=$ $-$$1/5$ has the lowest energy.
The state with $q_a$$=$ $-$$1/4$
is higher in energy by about 0.2 meV/Mn.
This behavior is similar to calculations without SO coupling (see Fig.~\ref{fig.Eq}).
Then, the AF2 state with $q_a$$=$ $-$$1/5$ appears to be nearly degenerate
with the AF1 state
(the energy difference is about $0.07$ meV/Mn, but the AF2 state is still low in energy).
Although these numerical values are probably on the verge of accuracy of our model analysis,
this tendency clearly shows that the SO interaction additionally stabilizes
AF1 state relative to the AF2 one. This seems to be reasonable: the conflict of
DM and isotropic exchange interactions, acting in the
opposite directions in one of the Mn sublattices, will
penalize the energy of the AF2 state. On the other hand, in the AF1 state, there is no such conflict
and the DM interactions will additionally minimize the energy of the $\uparrow \uparrow \downarrow \downarrow$ structure,
obtained without SO interaction. Apparently, the main reason why the energy of the AF2 phase is still
slightly lower than that of the AF1 phase is related to the fact that the isotropic exchange interactions
are overestimated in our approach (see Fig.~\ref{fig.J}) and, therefore, the effect of the SO coupling and
DM interactions on the magnetic structure is underestimated.

  Finally, we evaluate the value of the FE polarization in the AF2 phase, using the
Berry phase formalism,\cite{EPolarization} which was adopted for
the low-energy model.\cite{antispiral}
The vector of polarization is
parallel to the ${\bf b}$ axis, in agreement with the experiment.\cite{Taniguchi}
Then, the value of $P_b$ can be estimated as $2.0$, $3.8$, and $4.4$ $\mu$C/m$^2$ for
$q_a$$=$ $-$$1/3$, $-$$1/4$, and $-$$1/5$, respectively. Thus, $P_b$ is small, but it is
well consistent with the fact that the DM interactions are also small and
can be regarded as a small perturbation to the
homogeneous spin-spiral texture, formed by isotropic exchange interactions.
Nevertheless, we would like to emphasize that the
DM interactions are essential for deforming the spin spiral, breaking the inversion
symmetry, and producing the finite value of $P_b$.

  The experimental polarization is an order of magnitude larger: $P_b \sim  50 \mu$C/m$^2$.\cite{Taniguchi}
However, it can be still regarded as a `small value' in comparison with many other
multiferroic systems.\cite{MF_review} Therefore, it is quite consistent with our main idea that
$P_b$ is a result of small perturbation of the magnetic structure, caused by the relativistic effects.

  There may be several reason why the
experimental value of $P_b$ is larger than the theoretical one:

  (i) The experimental $P_b$ may also include
some lattice effects, in response to lowering of the magnetic symmetry;\cite{Shanavas}

  (ii) A unusual aspect of MnWO$_4$ is that the inversion symmetry is broken in the
``high-temperature'' phase, while the low-temperature phase remains centrosymmetric
(for comparison, the situation in perovskite manganites is exactly the opposite).\cite{MF_review}
Thus, in addition to the FE polarization, one should find some mechanism, which would stabilize
the phase AF2 over the phase AF1 at finite $T$. This mechanism may involve some temperature effects.

  (iii) The low-energy model is designed for the semi-quantitative analysis.
It is not always possible to expect a good quantitative agreement with the experimental data, because some
ingredients can be missing in the model. In principle, the GGA$+$$U$ calculations are also of a
semi-quantitative level, because the value of $P_b$ depends on the adjustable parameter $U$.\cite{remark1}
In our model analysis, the small value of $P_b$ may be again related to the overestimation of
isotropic exchange interactions (see Fig.~\ref{fig.J}). Then, the deformation of the spin-spiral
texture, caused by the relativistic effects, is underestimated. Therefore, $P_b$ is also
underestimated. From this point of view, the behavior of isotropic exchange interactions,
the total energies, and the FE polarization is consistent with each other,
and the main efforts towards quantitative description of the FE polarization and the phase
diagram of MnWO$_4$ should be concentrated on the correct description of isotropic exchange interactions.

\section{\label{sec:summary} Summary}

  In summary, we have provided the microscopic explanation for the origin of
FE activity in MnWO$_4$.
The multiferroicity in this compound is caused by the conflict of DM and
isotropic exchange interactions in the AF2 state.
Thus, MnWO$_4$ is multiferroic not simply because of the spiral structure.
It is essential to have conflicting interactions, which deform the spin spiral,
break the inversion symmetry and, thus, give rise to the FE activity.

  \textit{Acknowledgements}.
This work is partly supported by the grant of the Ministry of
Education and Science of Russia N 14.A18.21.0889.


\section*{Appendix I - Flavor associated energies for a two level system}

The time evolution of a quantum system of well-defined flavor quantum numbers described by the state vectors $\nu^{e}$ and $\nu^{\mu}$ respectively related to electron and muon neutrinos is given by
\begin{equation}
\left(\begin{array}{l}\nu^e_{(t)} \\ \nu^{\mu}_{(t)} \end{array}\right) =
U \, \left(\begin{array}{ll} e^{-i\, E_{\mbox{\tiny $1$}} t} & 0 \\ 0  &  e^{-i\, E_{\mbox{\tiny $2$}} t} \end{array}\right)\,
\left(\begin{array}{l}\nu_1 \\ \nu_2 \end{array}\right) =
U \, \left(\begin{array}{ll} e^{-i\, E_{\mbox{\tiny $1$}} t} & 0 \\ 0  &  e^{-i\, E_{\mbox{\tiny $2$}} t} \end{array}\right)
\, U^{\dagger} \left(\begin{array}{l} \nu^e_{(0)} \\ \nu^{\mu}_{(0)}  \end{array}\right),
\label{101A}
\end{equation}
where $\nu_{\mbox{\tiny $1$}}$ and $\nu_{\mbox{\tiny $2$}}$ are the mass eigenstates with well-defined energies, $E_{s} = \sqrt{p^{\mbox{\tiny $2$}} + m^{\mbox{\tiny $2$}}_{s}}$, with $s = 1,\, 2$, and the matrix $U$ parameterizes the mixing relation as
\begin{equation}
\left(\begin{array}{l}\nu^e_{(0)} \\ \nu^{\mu}_{(0)} \end{array}\right) =
U \, \left(\begin{array}{l}\nu_1 \\ \nu_2 \end{array}\right) =
\left(\begin{array}{rr} \cos\bb{\theta} & \sin\bb{\theta} \\ -\sin\bb{\theta}  & \cos\bb{\theta} \end{array}\right)
\left(\begin{array}{l}\nu_1 \\ \nu_2 \end{array}\right),
\label{req00A}
\end{equation}
where $\theta$ is the mixing angle.
Since the Hamiltonian of the system in the mass eigenstate basis can be extracted from Eq.~(\ref{req00A}) as $H = Diag\{E_{\mbox{\tiny $1$}},\, E_{\mbox{\tiny $2$}}\}$, the flavor projection operators can be easily defined as
\begin{equation}
M^{e}_{(t)} = |\nu^{e}_{(t)} \rangle \langle \nu^{e}_{(t)}| =
\left[
\begin{array}{ll}
\cos^{\mbox{\tiny $2$}}\bb{\theta} & \sin\bb{\theta} \,\cos\bb{\theta} \,e^{-i\,\Delta E\, t}\\
\sin\bb{\theta}\, \cos\bb{\theta} \, e^{i\,\Delta E\, t} &  \sin^{\mbox{\tiny $2$}}\bb{\theta}
\end{array}
\right]
\label{req01A}
\end{equation}
and
\begin{equation}
M^{\mu}_{(t)} = |\nu^{\mu}_{(t)} \rangle \langle \nu^{\mu}_{(t)}|  =
\left[
\begin{array}{ll}
\sin^{\mbox{\tiny $2$}}\bb{\theta} & -\sin\bb{\theta} \,\cos\bb{\theta} \,e^{-i\,\Delta E\, t}\\
-\sin\bb{\theta}\, \cos\bb{\theta} \, e^{i\,\Delta E\, t} &  \cos^{\mbox{\tiny $2$}}\bb{\theta} \end{array}
\right]
\label{req01B}
\end{equation}
where $\Delta E = E_{\mbox{\tiny $1$}} - E_{\mbox{\tiny $2$}}$ and it can be verified that $M^{e}_{(t)} + M^{\mu}_{(t)}  = \mathbf{1}$.

Thus the temporal evolution of a flavor eigenstate can be described by
\begin{equation}
|\nu^{e,\mu}_{(t)} \rangle = (M^{e}_{(0)} + M^{\mu}_{(0)}) |\nu^{e,\mu}_{(t)} \rangle =
 \langle \nu^{e}_{(0)} |\nu^{e,\mu}_{(t)} \rangle \, |\nu^{e}_{(0)} \rangle
+
\langle \nu^{\mu}_{(0)} |\nu^{e,\mu}_{(t)} \rangle \, |\nu^{\mu}_{(0)} \rangle,
\label{req01C}
\end{equation}
and the supposedly relevant measurable quantities, or {\em observables}, of the closed quantum system can be summarized by the the flavor-{\em averaged} energies,
\begin{eqnarray}
E^{e,\mu}_{(t)} &=&  \langle \nu^{e,\mu}_{(t)} |H| \nu^{e,\mu}_{(t)}\rangle,
\label{req02}
\end{eqnarray}
that result in time-independent quantities,
\begin{eqnarray}
E^{e}_{(t)} &=& E^{e}_{(0)} = \bar{E} + (1/2) \Delta E\, \cos\bb{2\theta},\nonumber\\
E^{\mu}_{(t)} &=& E^{\mu}_{(0)} = \bar{E} - (1/2) \Delta E\, \cos\bb{2\theta},
\label{req02A}
\end{eqnarray}
with $\bar{E} = (1/2)(E_{\mbox{\tiny $1$}} + E_{\mbox{\tiny $2$}})$, and by the time-oscillating flavor probabilities,
\begin{eqnarray}
\mathcal{P}_{\alpha\rightarrow\beta} \bb{t} &=& Tr\{M^{\beta}_{(0)}\, M^{\alpha}_{(t)}\}, ~~~~\alpha,\,\beta = e,\, \mu,
\label{req02BB}
\end{eqnarray}
that result in
\begin{eqnarray}
\mathcal{P}_{e\rightarrow e} \bb{t} &=& P_{\mu\rightarrow\mu} \bb{t} = |\langle \nu^{e}_{(\mbox{\tiny $0$})}|\nu^{e}_{(t)}\rangle|^{\mbox{\tiny $2$}} = |\langle \nu^{\mu}_{(\mbox{\tiny $0$})}|\nu^{\mu}_{(t)}\rangle|^{\mbox{\tiny $2$}}  = 1 - \sin^{\mbox{\tiny $2$}}\bb{2\theta} \,\sin^{\mbox{\tiny $2$}}\left(\frac{\Delta E}{2}\, t\right),
\label{req02B}
\end{eqnarray}
and
\begin{eqnarray}
\mathcal{P}_{e\rightarrow \mu}\bb{t} &=& P_{\mu \rightarrow e} \bb{t}= |\langle \nu^{\mu}_{(\mbox{\tiny $0$})}|\nu^{e}_{(t)}\rangle|^{\mbox{\tiny $2$}}  = |\langle \nu^{e}_{(\mbox{\tiny $0$})}|\nu^{\mu}_{(t)}\rangle|^{\mbox{\tiny $2$}} = \sin^{\mbox{\tiny $2$}}\bb{2\theta} \,\sin^{\mbox{\tiny $2$}}\left(\frac{\Delta E}{2}\, t\right),
\label{req02C}
\end{eqnarray}
that are interpreted as the probabilities of $e(\mu)$-flavor states produced at time $t_{\mbox{\tiny $0$}}$ be measured as $e(\mu)$-flavor states or be converted into $\mu (e)$-flavor states after a time interval $t - t_{\mbox{\tiny $0$}} \sim t - 0 \sim t$.

Now let us suppose that the density matrix of a composite quantum system of two neutrino flavor states is given by
\begin{equation}
\hat{\rho}\bb{t} \equiv \hat{\rho} = w_e \, M^{e}_{(t)} + w_{\mu} \, M^{\mu}_{(t)},
\label{req05}
\end{equation}
with $w_{e} + w_{\mu} = 1$.
One easily finds that the re-defined probabilities of measuring the electron and muon flavor eigenstates at time $t$ are given by
\begin{eqnarray}
P^{e}_{(t)}   = \mbox{Tr}\{M^{e}_{(0)}\, \hat{\rho}\} &=& w_e \, \mbox{Tr}\{M^{e}_{(0)}\,M^{e}_{(t)}\} + w_{\mu} \, \mbox{Tr}\{M^{e}_{(0)}\,M^{\mu}_{(t)}\} = \nonumber\\
                        &=& w_e \, \mathcal{P}_{e\rightarrow e}\bb{t} + w_{\mu} \, \mathcal{P}_{\mu\rightarrow e}\bb{t},\\
P^{\mu}_{(t)} = \mbox{Tr}\{M^{\mu}_{(0)}\, \hat{\rho}\} &=& w_e \, \mbox{Tr}\{M^{\mu}_{(0)}\,M^{e}_{(t)}\} + w_{\mu} \, \mbox{Tr}\{M^{\mu}_{(0)}\,M^{\mu}_{(t)}\} = \nonumber\\
                        &=& w_e \, \mathcal{P}_{e\rightarrow \mu}\bb{t} + w_{\mu} \, \mathcal{P}_{\mu\rightarrow \mu}\bb{t},
\label{req06}
\end{eqnarray}
where we have used the results from Eqs.~(\ref{req02B}-\ref{req02C}).
One also easily notices that
\begin{eqnarray}
P^{e}_{(t)} +  P^{\mu}_{(t)} &=&
w_{e} (\mathcal{P}_{e\rightarrow e}\bb{t} + \mathcal{P}_{e\rightarrow \mu}\bb{t})
+
w_\mu (\mathcal{P}_{\mu\rightarrow e}\bb{t} + \mathcal{P}_{\mu\rightarrow \mu}\bb{t})
= w_{e} + w_\mu = 1
\label{req06B}
\end{eqnarray}
and that the properties of a statistical mixture are immediate.
It leads to a reinterpretation of the energy related to each flavor quantum number.

The standard {\em total} averaged energy for a composite quantum system is defined through the density matrix as
\begin{eqnarray}
E_{(t)} = \mbox{Tr}\{H \, \hat{\rho}\} &=& w_{e} \, \mbox{Tr}\{H \, M^{e}_{(t)}\} + w_{\mu} \, \mbox{Tr}\{H \, M^{\mu}_{(t)}\}\nonumber\\
                            &=&   w_{e} \, E^{e}_{(t)} + w_{\mu} \, E^{\mu}_{(t)},
\label{req07}
\end{eqnarray}
from which one can notice the explicit dependence on the flavor-{\em averaged} energies, $E^{e,\mu}_{(t)}$, recovered from Eq.~(\ref{req02}).
In this context $E^{e}_{(t)}$ and $E^{\mu}_{(t)}$ are respectively decoupled from the statistical weights $w_{\mu}$ and $w_{e}$.
It just ratifies our previous arguments that such flavor energies are noway correlated with the flavor probabilities from Eq~(\ref{req07}), $P^{e}_{(t)}$ and $P^{\mu}_{(t)}$ since both of them depend simultaneously on both statistical weights, $w_{\mu}$ and $w_{e}$.
Thus the arguments that assert the ambiguity and the insufficiency in defining the flavor eigenstate averaged energies through $E^{e,\mu}_{(t)}$ are maintained.

After simple mathematical manipulations involving the definitions from Eq.~(\ref{req01B}) and the probabilities from Eq.~(\ref{req06}), one easily finds that
\begin{eqnarray}
M^{\mu}_{(0)} \, \hat{\rho}\, M^{\mu}_{(0)} &=& (w_e \, \mathcal{P}_{e\rightarrow \mu} \bb{t}+ w_{\mu} \, \mathcal{P}_{\mu\rightarrow \mu}\bb{t}) M^{\mu}_{(0)} = P^{\mu}_{(t)}\, M^{\mu}_{(0)},\nonumber\\
M^{e}_{(0)} \, \hat{\rho}\, M^{e}_{(0)}     &=& (w_e \, \mathcal{P}_{e\rightarrow e}\bb{t} + w_{\mu} \, \mathcal{P}_{\mu\rightarrow e}\bb{t}) M^{e}_{(0)} = P^{e}_{(t)}\, M^{e}_{(0)},
\label{req08}
\end{eqnarray}
Observing the cyclic properties of the trace, the flavor-{\em weighted} energies can be defined as
\begin{eqnarray}
\epsilon^{e,\mu}_{(t)} = \mbox{Tr}\{M^{e,\mu}_{(0)}\,H\, M^{e,\mu}_{(0)} \, \hat{\rho} \} &=& \mbox{Tr}\{H\, M^{e,\mu}_{(0)} \, \hat{\rho}\, M^{e,\mu}_{(0)}\}
= P^{e,\mu}_{(t)}\mbox{Tr}\,\{M^{e,\mu}_{(0)}\,H\} =
P^{e,\mu}_{(t)}\, E^{e,\mu}_{(0)},
\label{req09}
\end{eqnarray}
which can be promptly compared with the previous definition through the relation
\begin{equation}
\frac{|\epsilon^{e,\mu}_{(t)} - w_{e,\mu} E^{e,\mu}_{(0)}|}{E^{e,\mu}_{(0)}} = |w_{e} - w_{\mu}| \sin^{\mbox{\tiny $2$}}\bb{2\theta}
\, \sin^{\mbox{\tiny $2$}}\left(\frac{\Delta E}{2}\,t\right).
\label{req10A}
\end{equation}

As described above, flavor-{\em weighted} energies establish a unique correspondence between flavor eigenstate energies and the statistical definitions of probabilities, $P^{e,\mu}_{(t)}$.
It can not be identified through the definition of flavor-{\em averaged} energies.

\section*{Appendix II - The von-Neumann entropy and quantum measurements}

The connection to the quantum measurement theory arises from assuming the von-Neumann entropy as an important entropy functional defined in terms of the density matrix by
\begin{equation}
S\bb{\hat{\rho}} = - \mbox{Tr}\{\hat{\rho}\bb{t}\, \ln\bb{\hat{\rho}\bb{t}}\},
\label{eq26Z}
\end{equation}
where we have set the multiplicative Boltzmann constant, $k_{B}$, equal to unity.
The entropy $S\bb{\hat{\rho}}$ quantifies the departure of a composite quantum system from a pure state by measuring its time-evolved degree of mixture.
The quantum measurements can also induce modifications on the the von-Neumann entropy of the system.
The entropy changes due to a {\em non-selective} measurement scheme described by {\em operations} parameterized by the projection operators $M^{\alpha}_{(0)}$ are given by \cite{Breuer}
\begin{equation}
\Delta S = S\bb{\hat{\rho}^{\prime}} - S\bb{\hat{\rho}} \geq 0,
\label{eq26A}
\end{equation}
where $\alpha$ is the relevant quantum number and $\hat{\rho}^{\prime} = \sum_{\alpha} {P^{\alpha}_{(t)} \hat{\rho}_{\alpha}}$ so that
\begin{equation}
S\bb{\hat{\rho}^{\prime}} =
S\left(\sum_{\alpha} {P^{\alpha}_{(t)} \hat{\rho}_{\alpha}}\right).
\label{eq26B}
\end{equation}
Since $\Delta S \geq 0$, the {\em non-selective} quantum measurement never decreases the von-Neumann entropy.
To quantify the relation between {\em selective} and {\em non-selective} levels of measurement, the mixing entropy described by
\begin{equation}
\delta S =
S\left(\sum_{\alpha} {P^{\alpha}_{(t)} \hat{\rho}_{\alpha}}\right) - \sum_{\alpha} {P^{\alpha}_{(t)} S\left(\hat{\rho}_{\alpha}\right)}.
\label{eq26C}
\end{equation}
gives the difference between the entropy of a system projected by a {\em non-selective} quantum measurement, $S\left(\sum_{\alpha} {P^{\alpha}_{(t)} \hat{\rho}_{\alpha}}\right)$, and the average of the entropies of the sub-ensembles $\hat{\rho}_{\alpha}$, described by $M^{\alpha}_{(0)}$.
For the {\em selective} measurement scheme, with $M^{\alpha}_{(0)}$ denoting the creation of a single-flavor (pure) state, the mixing entropy is reduced to
\begin{equation}
\delta S = S\left(\sum_{\alpha} {P^{\alpha}_{(t)} \hat{\rho}_{\alpha}}\right),
\label{eq26D}
\end{equation}
since $S\left(\hat{\rho}_{\alpha}\right) = 0$.
The {\em non-selective} measurement described by observables like flavor-{\em weighted} energies always modifies the von-Neumann entropy by $\Delta S$ (c. f. Eq.~(\ref{eq26A})).
In the particular case of an ensemble described as a pure state, the {\em selective} measurement of flavor-{\em averaged} energies which are not expressed in terms of flavor conversion probabilities does not modify the von-Neumann entropy.
From the theoretical perspective, the von-Neumann entropy is therefore an auxiliary variable in distinguishing the measurement procedure and in classifying the measurement interventions as {\em selective} and {\em non-selective} ones, which however, as noticed above, are not complementary concepts.

All the above defined entropies satisfy some set of inequalities \cite{Breuer} which have been extensively used in different forms in the framework of quantum information theory and quantum entanglement.