\section{Introduction}\label{sec:section1}
Since the accelerating expansion of the universe was revealed by observations of type Ia supernovae (SN) ~\citep{Riess:1998, Perlmutter:1999}, SN plays an essential role to study the expansion history of the Universe and the nature of dark energy. However, due to the limited intrinsic luminosity and the extinction from the interstellar medium, the maximum redshift of the SN we can currently observe is about $1.7$. Approximately, for most other popular probes at low redshift, such as baryon acoustic oscillations (BAO) and weak lensing, the bottleneck of redshift range in the near future observations is about $z\sim2$. Whereas fluctuations of the cosmic microwave background (CMB) provide cosmological information at the last scattering surface ($z\sim1089$). Therefore, we do not have many effective methods to observe the evolution of our Universe at $2<z<1089$, usually referred as the ``cosmological desert''.

Gamma-ray bursts (GRB) are the most energetic explosions after the big bang in the Universe. The equivalent isotropic energy radiated in a few seconds can reach up to $10^{48}\sim10^{53}$ erg (see e.g.,~\citep{Piran:1999, Meszaros:2002, Meszaros:2006, Kumar and Zhang:2015} for recent reviews). Thanks to the extreme brightness and the immunity to dust extinction of high-energy photons in the gamma-ray band, GRB are detectable up to redshift $z\sim 10$ ~\citep{Salvaterra:2009, Tanvir:2009, Cucchiara:2011, Tanvir:2013}. Therefore, in the literatures, GRB are widely proposed, as a complementary data of SN, to trace the Hubble diagram of the Universe at high redshifts. However, different from the consistent luminosities of SN which allows us to qualify SN as an ideal distance indicator, the central engine mechanism of explosions of GRB has not yet clearly been understood. Therefore, drawing cosmological implications from GRB observations is quite intractable. Recently, many methods have been proposed to achieve some progresses, in order to treat GRB as the distance indicator ~\citep{Schaefer:2002, Dai:2004tq, Ghirlanda:2004a, Liang:2005, Firmani:2005, Schaefer:2007, Liang:2008, Li:2008, Wei:2009, Wei:2010, Liu:2015, Li:2015, Lin:2016,Wang:2015}.

Luminosity or energy correlations of GRB are empirical connections between measurable properties of the prompt gamma-ray emission and the luminosity or energy of GRB. In recent years, several empirical luminosity correlations have been statistically proposed from observations, such as various popular two-variable relations: the correlation between spectrum lag and isotropic peak luminosity ($\tau_{\rm lag}-L$)~\citep{Norris:2000}; the correlation between time variability and isotropic peak luminosity ($V-L$)~\citep{Fenimore:2000}; the correlation between the peak energy of $\nu F_{\nu}$ spectrum and isotropic equivalent energy ($E_{\rm p}-E_{\rm iso}$, the Amati relation)~\citep{Amati:2002}; the correlation between peak energy and collimation-corrected energy ($E_{\rm p}-E_{\gamma}$)~\citep{Ghirlanda:2004b}; the correlation between peak energy and isotropic peak luminosity ($E_{\rm p}-L$)~\citep{Schaefer:2002,Yonetoku:2004}; the correlation between minimum rise time of light curve and isotropic peak luminosity ($\tau_{\rm RT}-L$)~\citep{Schaefer:2007}; and some well-known multi-variable relations, such as the correlation between $E_{\rm iso}$, $E_{\rm p}$, and the break time of the optical afterglow light curves $t_{\rm b}$~\citep{Liang:2005}. The general expression for luminosity correlations is
\begin{equation}\label{eq1}
  y=a+bx~,
\end{equation}
where $x=\log E_{\rm p}$ (or $\log V$, $\log\tau_{\rm lag}$, $\log\tau_{\rm RT}$), $y=\log L$ (or $\log E_{\rm iso}$, $\log E_{\gamma}$). Determinations of coefficients $a$ and $b$ of correlations could tell us whether these empirical luminosity correlations are independent on the redshift, which is very crucial for treating GRB as the distance indicator to investigate the evolution of our Universe. In addition, some new models for GRB have been proposed in resent years, such as the spectro-temporal multi-component model \citep{Guiriec:2011,Guiriec:2013,Guiriec:2015a,Guiriec:2015b,Guiriec:2016a,Guiriec:2016b}. Interestingly, this model could result in a new time-resolved relation: the relation between the luminosity of the non-thermal component, $L^{\rm nTH}_{\rm i}$, and its corresponding $\nu F_{\nu}$ spectral peak energy in the rest frame, $E^{\rm NT,rest}_{\rm peak,i}$ (i.e., $L^{\rm nTH}_{\rm i}-E^{\rm NT,rest}_{\rm peak,i}$ relation).

More importantly, determinations of coefficients can also show the reliability of the calibration of GRB luminosity correlations directly. In the literatures, there are two methods to calibrate the luminosity correlations: ``self-calibration'' method, which is based on the global fitting technique using all GRB data to constrain the cosmological parameters and the coefficients of luminosity correlations simultaneously~\citep{Schaefer:2002, Li:2008}; ``low-redshift calibration'' method, which is using the low-redshift GRB data to obtain the constraints on coefficients of luminosity correlations, and directly extrapolate them to high redshifts to constrain the cosmological parameters using the high-redshift GRB data~\citep{Liang:2008, Wei:2009, Wei:2010, Ding:2015, Liu:2015, Lin:2016, Wang:2016}.

Apparently, these two methods are seriously dependent on the assumption that empirical luminosity correlations are universal and do not evolve with respect to redshift. In practice, for the Amati relation, \citet{Li:2007} first revealed that the systematically significant variations of the intercept $a$ and the slope $b$ with the cosmological redshift. Later, \citet{Basilakos:2008} investigated the same issue for several other luminosity relations. Due to the limited quality of GRB samples, there was no statistically significant evidence for evolution of the calibration parameters. In addition, as a cosmological probe, the evolution of GRB luminosity correlations may result in overestimate or underestimate of cosmological parameters~\citep{Dainotti:2013}. More recently, in~\citet{Lin:2015}, they re-investigated the Amati relation using low-redshift ($z < 1.4$) and high-redshift ($z > 1.4$) GRB, and found that the coefficients of Amati relation in low-redshift GRB differs from those of high-redshift GRB at more than $3\sigma$ confidence level. Therefore, they finally concluded that long gamma-ray bursts might not be standard candles. Following on this topic, \citet{Lin:2016} further examined the possible redshift dependence of several other luminosity correlations and found that only for the $E_{\rm p} - E_\gamma$ relation, the low-redshift GRB could give the similar constraints on the coefficients with those from the high-redshift GRB within $1\sigma$ confidence level.

In this paper, we choose the Amati relation as an example to investigate the evolutions and calibrations of GRB luminosity correlation with the latest GRB observational data~\citep{Liu:2015}. Firstly, we divide the whole GRB sample into two bins ($z<1.4$ and $z>1.4$). In each bin, we use the GRB data to constrain the coefficients in the $\Lambda$CDM framework to check whether they are independent on the redshift. Then, we propose an extended Amati relation by introducing two extra terms to characterize the evolutions with respect to the redshift. Secondly, we carefully check the reliability of the calibration method. In practice, we use the low-redshift GRB data ($z<1.4$) to calibrate the coefficients of GRB luminosity correlation and directly extrapolate them to high-redshift events. Then we use this ``calibrated'' Amati relation to constrain the matter density parameter $\Omega_{\rm m}$ with the high-redshift GRB data. Finally, we simulate several samples of GRB data with high precisions to investigate that in which conditions the GRB data calibrated by the low-redshift calibration method can be used for cosmological studies.

This paper is organized as follows: Section \ref{sec:evolution} and Section \ref{sec:calibration} are dedicated to present the method and numerical results on the analysis of evolution and calibration of GRB luminosity correlation. Finally, conclusions and discussions are presented in section \ref{sec:summary}.

\section{Evolutions of Luminosity Correlations}\label{sec:evolution}

\subsection{Amati Relation \& GRB Data}\label{sec:section2.1}

Although for recent years, many empirical luminosity correlations have been statistically concluded from long GRB observations, the Amati correlation is the most widely used one among them for cosmological studies. The original version of the Amati relation is expressed as
\begin{align}
\log{\frac{E_{\text{iso}}}{\text{erg}}}&=a+b\log{\frac{E_{{\rm p},i}}{300~\text{keV}}}~,\label{equ:relation}
\end{align}
where $E_{{\rm p},i}=E_{\rm p}(1+z)$ to correct the redshift dilation of the spectrum~\citep{Wang:2011}. The isotropic equivalent energy $E_{\rm iso}$ can be calculated from the bolometric fluences $S_{\rm bolo}$~\citep{Schaefer:2007}:
\begin{equation}\label{equ:Eiso}
E_{\rm iso}=4\pi d^{2}_{L}S_{\rm bolo}\left(1+z \right)^{-1}~.
\end{equation}
The uncertainty of $E_{\rm iso}$ propagates from the uncertainties of $S_{\rm bolo}$ and $d_{\rm L}$. $S_{\rm bolo}$ is calculated from the observed peak photon flux in the rest frame $1-10,000$ keV energy band by assuming the Band spectrum~\citep{Band:1993}. To calculate the luminosity distance $d_{\rm L}$, we assume a flat concordance $\Lambda$CDM model with the Hubble constant: $H_0=67.8~\rm km~s^{-1}Mpc$ and the matter density parameter: $\Omega_{\rm m}=0.308$ obtained from the latest $Planck$ 2015 results~\citep{Ade:2015}. Hence, $d_{\rm L}$ in Eq.~(\ref{equ:Eiso}) can be written as
\begin{equation}\label{equ:luminosity}
d_{\rm L}(z)=\frac{c(1+z)}{H_0}\int^{z}_{0}\frac{dz}{\sqrt{\Omega_{\rm m}(1+z)^3+(1-\Omega_{\rm m})}}~.
\end{equation}

\begin{table*}
\centering
\caption{$1\,\sigma$ Constraints on the intrinsic scatters ($\sigma_{\rm int}$) and luminosity correlation parameters ($a$, $b$, $\alpha$ and $\beta$) in the flat $\Lambda$CDM framework for the original Amati relation and extended Amati relation respectively. For comparison, we also show the constraints using the global fitting method.}\label{tab:parameters_LCDM}
\begin{tabular}{llccccccc}
\hline\hline
relation   & bins & $\sigma_{\rm int}$ & N & $\chi^2_{\rm min}$ & $a$ & $b$ & $\alpha$ & $\beta$\\
\hline
&low-$z$& $0.40\pm0.05$ & $59$ & $55.51$ & $52.75\pm0.06$ & $1.60\pm0.10$ & &\\
Amati relation:&high-$z$& $0.32\pm0.04$ & $79$ & $76.18$ & $52.95\pm0.05$ & $1.30\pm0.12$ & &\\
&Global Fit& $0.35\pm0.03$ & $138$ & $134.83$ & $52.62\pm0.07$ & $1.49\pm0.07$ & &\\
\hline
Extended &low-$z$& $0.40\pm0.05$ & $59$ & $55.15$ & $52.42\pm0.35$ & $1.46\pm0.29$ & $0.72\pm0.75$ & $0.14\pm0.64$\\
Amati relation: &high-$z$& $0.32\pm0.04$ & $79$ & $74.06$ & $53.42\pm0.46$ & $0.08\pm1.01$ & $-0.67\pm0.64$ & $1.73\pm1.42$\\
&Global Fit& $0.35\pm 0.03$ & $138$ & $134.65$ & $52.52\pm0.16$ & $1.59\pm0.18$ & $0.33\pm 0.40$ & $-0.31\pm0.32$\\
\hline\hline
\end{tabular}
\end{table*}

\begin{figure}[t]
	\centering
	\includegraphics[width=0.45\textwidth]{Amati_relation.eps}
	\caption{The original Amati relation for all GRB data. The blue points with error bars are GRB data. The red line denotes the theoretical prediction from the best-fit values $a=52.82$ and $b=1.56$, which are derived from maximizing the D'Agostini's likelihood.}\label{fig:Amati_relation}
\end{figure}

To test the possible redshift dependence of the Amati luminosity correlation, we use the GRB sample compiled in~\citet{Liu:2015}, which includes 138 well-measured GRB in the redshift range: $z\in [0.0331,~8.1]$. In general, we examine the evolution of the coefficients $a$ and $b$ according the following steps: Firstly, we get the equivalent isotropic energy $E_{\rm iso}$ for each GRB at redshift $z$ with the distance in Eq.~(\ref{equ:luminosity}) taken into account. Then, we divide all observed 138 GRB sample into two redshift bins: [0.0331, 1.4], [1.4, 8.1]. The reason of choosing $z = 1.4$ as the threshold is that the maximal redshift of SNe Ia data samples we are using for cosmological studies is usually around $z=1.4$. The Universe below $z = 1.4$ has already been well studied by using the SNe Ia datasets \citep{Amanullah:2010, Suzuki:2012, Betoule:2014}. Therefore, we treat the GRB data with $z < 1.4$ and $z > 1.4$ as the low-z and high-z samples, respectively. In these two redshift bins, there are 59 and 79 GRB data, respectively. And then, we constrain the luminosity correlation coefficients $a$ and $b$ by maximizing the D'Agostini's likelihood~\citep{DAgostini:2005}:
\begin{align}\label{D_like}
\nonumber \mathcal{L}_D(a,b,\sigma_{\rm int})\propto
&\prod_i\frac{1}{\sqrt{\sigma_{\rm int}^2+\sigma_{y_i}^2+b^2\sigma_{x_i}^2}}\\
&\times \exp\left[-\frac{(y_i-a-bx_i)^2}{2(\sigma_{\rm int}^2+\sigma_{y_i}^2+b^2\sigma_{x_i}^2)}\right],
\end{align}
or, equivalently, by minimizing the $\chi^2$:
\begin{align}
\nonumber\chi^{2}_{D}(a,b,\sigma_{\rm int})=&\sum_{i}\ln \left(\sigma^{2}_{\rm int}+\sigma^{2}_{y_{i}}+b^{2}\sigma_{x_{i}}^{2}\right)\\
&+\sum_{i}\frac{\left(y_{i}-a-bx_{i}\right)^{2}}{\sigma^{2}_{\rm int}+\sigma^{2}_{y_{i}}+b^{2}\sigma_{x_{i}}^{2}},
\end{align}
where $\sigma_{\rm int}$ is the intrinsic scatter which represents any other unknown uncertainties except for the observational statistical ones. Here, in order to test the feasibility of D'Agostini's likelihood, we constrain the two parameters of Amati relation use all GRB data (138 GRB) in the flat $\Lambda$CDM framework firstly. The minimum chi-squared value is $\chi^2_{\rm min}=134.91$ and $\chi^2_{\rm min}/{\rm d.o.f.}=0.98$. We also plot the comparison between the GRB data and the theoretical prediction of the best-fit values in Fig. \ref{fig:Amati_relation}. All these results show that the D'Agostini’s likelihood is feasible to constrain parameters. Finally, after getting constraints on $a$ and $b$, in Fig.\ref{fig:testnoevo} we plot the 2 dimensional marginalized constraints in the plane of ($a$, $b$) from the low-redshift (blue contours) and high-redshift (red contours) GRB data, respectively.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.43\textwidth]{contour_LCDM5_2-5_0.eps}
	\caption{The 2-dimensional marginalized distributions with 1$\sigma$, 2$\sigma$, and 3$\sigma$ contours for the coefficients $a$, $b$, $\alpha$, and $\beta$ for the original Amati relation (the upper right panel) and the extended Amati relation (other six panels), respectively. The blue and red contours denote the results on the parameters of the GRB relations from the low-redshift and high-redshift GRB subsamples, respectively. The central points represent the best-fit values.}\label{fig:testnoevo}
\end{figure}

\begin{figure*}[t]
	\centering
	\includegraphics[width=0.8\textwidth]{Omega_m_LCDM_glbfit5_2-5_0.eps}
	\caption{1-dimensional and 2-dimensional marginalized constraints on $\Omega_{\rm m}$ and $\alpha$ for the original Amati relation (the upper right panel) and the extended Amati relation (other three panels). For the extended Amati relation, blue lines of the lower left three panels represent the results from global fitting method. The dotted red lines in the left two panels represent $\alpha=0.72$ (see Table \ref{tab:parameters_LCDM}) which is the best-fit value obtained from the low-redshift GRB data in the low-$z$ calibration method. The constraint on $\Omega_{\rm m}$ from the low-redshift calibration method are shown in the lower right panel (the red dashed line, the green dash dot line, and the black dotted line). For comparison, in the upper right panel we also show the results obtained from the low-redshift calibration method and the global fitting method for the original Amati relation, where the blue line is the result of global fitting method, the red dashed line and green dash dot line are the result from low-redshift calibration method.}\label{fig:Omm}
\end{figure*}

In our analysis, we use {\it emcee}\footnote{https://pypi.python.org/pypi/emcee}~introduced by \citet{Foreman-Mackey:2012}, a Python module that uses the Markov chain Monte Carlo (MCMC) method to get the best-fit values and their uncertainties of parameters $a$, $b$ and $\sigma_{\rm int}$ by generating sample points of the probability distribution. The best-fit values with 1$\sigma$ errors are shown in Table \ref{tab:parameters_LCDM}. We also plot the results in Fig.~\ref{fig:testnoevo} (the upper right panel), in which the blue and red contours denote the results on the parameters of Amati relation from the low-redshift and high-redshift GRB subsamples, respectively. Apparently, these two constraints are not match very well. The best fit values are more than $3\sigma$ away. When considering the uncertainties, the tension between two fitting results is less than $2\sigma$ confidence level, which implies that the parameters of Amati relation might be redshift-dependent. The result is similar with some other works \citep{Lin:2015}.

\subsection{Redshift Evolution of Amati Relation}\label{sec:section2.2}

As the statistically obvious evolution of the luminosity correlation has been revealed, we should be careful to use the GRB data as the distance indicator for cosmological studies. Because all these investigations are based on the assumption that the luminosity correlation does not vary with respect to the redshift, i.e., both $a$ and $b$ in Eq.~(\ref{equ:relation}) should be constant. To solve this problem, inspired by the study for the evolution of light-curve fitting parameters in SN observations~\citep{Wang:2013, Wang:2014} and that for the evolution of the GRB intrinsic luminosity~\citep{Yu:2015}, we introduce two extra redshift-dependent terms to characterize the redshift evolution of the luminosity correlation. Since most redshifts of GRB are much greater than unity, rather than the linear expression in~\citet{Wang:2013, Wang:2014} and the logarithmic parametrization in~\citet{Yu:2015}, we select two mild formulas
\begin{equation}\label{equ2}
  a\rightarrow A=a+\alpha\frac{z}{1+z}~~,~~~b\rightarrow B=b+\beta\frac{z}{1+z}
\end{equation}
to avoid extreme results at high redshifts. Therefore, the extended Amati relation can be expressed as
\begin{align}
\log{\frac{E_{\text{iso}}}{\text{erg}}}&=\bigg(a+\frac{\alpha z}{1+z}\bigg)+\bigg(b+\frac{\beta z }{1+z}\bigg)\log{\frac{E_{\rm p,i}}{300~\text{keV}}},\label{equ:newrelation}
\end{align}
where $a$,$b$,$\alpha$ and $\beta$ are new four coefficients in the relation which should be obtained from the GRB data.

Following the same procedures for analyzing the evolutions of $a$ and $b$ presented in the subsection~\ref{sec:section2.1}, we also investigate the evolutions of luminosity correlation coefficients $a$, $b$, $\alpha$ and $\beta$ using the GRB sample in the same low-redshift and high-redshift subsamples, respectively. The best-fit values with 1$\sigma$ errors are shown in Table \ref{tab:parameters_LCDM}. Fig. \ref{fig:testnoevo} (the lower left six panels) shows the results regarding evolutions of luminosity correlation coefficients ($a$, $b$, $\alpha$ and $\beta$). Since we have two parameters to describe the redshift dependence, compared with the result of original Amati relation, the tension between constraints of parameters from low-$z$ and high-$z$ GRB data are alleviated, most of which are consistent within $1\sigma$ confidence level. Of course, these two extra free parameters could bring the large uncertainties and make the constraints weaker.

\section{Calibrations of Luminosity Correlations}\label{sec:calibration}

\subsection{Low Redshift Calibration}\label{sec:section3.1}

As shown in section \ref{sec:section2.2}, the tension between the constraints on parameters of the original Amati relation from the low-redshift and high-redshift GRB subsamples is more than $1\sigma$ confidence level, while the constraints are consistent within $1\sigma$ confidence level for the extended Amati relation. However, this does not mean that we can safely use the GRB data for cosmological studies. We need to check the calibration of GRB luminosity correlations further. Usually, we have two calibration methods. One is to calibrate luminosity correlations with low-redshift GRB samples in the context of the flat $\Lambda$CDM model and then extrapolate these obtained coefficients to high redshift range. The other method is the global fitting analysis using the MCMC technique in which luminosity correlation parameters and cosmological parameters are simultaneously fitted on the same weight in the context of the flat $\Lambda$CDM model \citep{Li:2008}. Therefore, here we use these two methods to calibrate the extended Amati relation and compare the constraints on the cosmological parameters.

In the low-redshift calibration method, we obtain the luminosity distance of the low-redshift GRB sample ($z < 1.4$) by using Eq.~(\ref{equ:luminosity}). And then we calibrate the luminosity correlation by maximizing the D'Agostini's likelihood (Eq.~\ref{D_like}) for the 59 low-redshift GRBs. The constraints on the luminosity correlation coefficients are shown in Tab. \ref{tab:parameters_LCDM}. Then we directly extrapolate these coefficients to high redshifts. Therefore, we construct the Hubble diagram for the 79 high-redshift GRBs based on the calibrated luminosity correlation and investigate cosmological implication from this Hubble diagram in the context of standard $\Lambda$CDM model. In this framework, the luminosity distance is determined by Eq.~(\ref{equ:luminosity}), and the corresponding distance modulus is
\begin{equation}
\mu(z)=5\lg\frac{d_L(z)}{\rm Mpc}+25.
\end{equation}
We use the numerical method to fit $\Omega_{\rm m}$ and the corresponding $\chi^2$ is
\begin{equation}\label{equ:SNe-glbfit-chi2}
\chi^2 = \sum_i\frac{[\mu_{\rm obs}(z_i)-\mu_{\rm th}(z_i;\Omega_{\rm m})]^2}{\sigma^2_{\mu,i}},
\end{equation}
where $\sigma_{\mu,i}$ is the total uncertainty of distance modulus for the $i$th GRB and it is propagated from the uncertainties of $S_{\rm bolo},~E_{\rm p,i},~E_{\rm iso}$
\begin{equation}\label{error}
\sigma^2_{\mu}=\left(\frac{5}{2}\right)^2\left[\sigma^2_{\log E_{\rm iso}}+\left(\frac{\sigma_{S_{\rm bolo}}}{S_{\rm bolo}\ln 10}\right)^2\right],
\end{equation}
where
\begin{align*}
&\sigma^2_{\log{E_{\rm iso}}}=~\sigma^2_a+\left(\frac{z\sigma_{\alpha}}{1+z}\right)^2+\left(\sigma_{b}\log\frac{E_{\rm p,i}}{300~\rm keV}\right)^2\\
&+\left(\frac{z\sigma_{\beta}}{1+z}\cdot\log\frac{E_{\rm p,i}}{300~\rm keV}\right)^{2}+
\left[\frac{\left(b+\frac{z\beta}{1+z}\right)\sigma_{E_{\rm p,i}}}{E_{\rm p,i}\ln{10}}\right]^2+\sigma^2_{\rm int}.
\end{align*}

Firstly, we consider the original Amati relation. In the upper right panel of Fig. \ref{fig:Omm} (red dashed line), we show the one-dimensional distribution of the matter density parameter $\Omega_{\rm m}$ obtained from the 79 high-redshift GRB data by using the low-redshift calibration method. If we use the best-fit values $a=52.75$ and $b=1.60$ to calibrate the high-redshift GRB data, we obtain the 68\% C.L. constraint on the matter density parameter:
\begin{equation}\label{eq:omm:low-z0}
  \Omega_{\rm m}=0.50\pm0.10~.
\end{equation}
Based on this result, we can see that this obtained constraint on $\Omega_{\rm m}$ is far away from the fiducial value of $\Omega_{\rm m}=0.308$, which is used to constrain the parameter $a$ and $b$ from the low-$z$ GRB data for calibration. The significance of this inconsistence is about $2\sigma$ confidence level.

Then, we use the extended Amati relation to do the calculations, following the same procedures. In the lower right panel of Fig. \ref{fig:Omm} (red dashed line), we also show the one-dimensional distribution of $\Omega_{\rm m}$ obtained from the 79 high-redshift GRB data by using the low-redshift calibration method. The obtained constraints on $\Omega_{\rm m}$ is:
\begin{equation}\label{eq:omm:low-z}
  \Omega_{\rm m}=0.24\pm0.12~,
\end{equation}
at 68\% confidence level, which is consistent with previous works \citep{Wang:2016}. This result is consistent with the fiducial value of the matter density parameter very well, which also implies that the extended Amati relation might be better than the original one.

\begin{figure}
	\centering
	\includegraphics[width=0.45\textwidth]{contour_LCDM_glbfit5_0.eps}
	\caption{The 2-dimensional marginalized distributions with 1$\sigma$, 2$\sigma$, and 3$\sigma$ contours for the coefficients $a$ and $b$ for the original Amati relation. The blue and green contours denote the results obtained from the low-$z$ GRB data in the low-redshift calibration method and from all GRB data in the global fitting method. The points represent the best-fit values of coefficients.}\label{fig:glbfit_LCDM}
\end{figure}

\subsection{Global Fitting Method}\label{sec:section3.2}

The second method is the global fitting method, in which we calculate the matter density parameter $\Omega_{\rm m}$ and luminosity correlation coefficients ($a$, $b$, $\alpha$ and $\beta$) simultaneously with the same weight. The constraints on these coefficients from all-redshift GRB data are also listed in Table \ref{tab:parameters_LCDM}.

Similar with the discussions in the section \ref{sec:evolution}, we start with the original Amati relation. In Table \ref{tab:parameters_LCDM} and Fig. \ref{fig:glbfit_LCDM}, we show that constraints on the parameters of Amati relation from the low-$z$ GRB data in the low-redshift calibration method and from all GRB data in the global fitting method are more or less consistent within $1\sigma$ confidence level. As we know, the constraint on the matter density parameter is strongly correlated with these coefficients of GRB relation. Therefore, the constraint on $\Omega_{\rm m}$ using the global fitting method is quite different from that using the low-redshift calibration method:
\begin{equation}
  \Omega_{\rm m}>0.45
\end{equation}
at $95\%$ confidence level, as shown in the blue line of the upper right panel of Fig. \ref{fig:Omm}. This large difference implies that the obtained constraint on $\Omega_{\rm m}$ using the low-redshift calibration method could be biased and overestimated. On the other hand, the constraint on $\Omega_{\rm m}$ using the global fitting method is far away from the current best-fit values from other observational data $\Omega_{\rm m}\sim0.3$, which might imply that the original Amati relation could not be the best relation to describe the GRB data.

Next, we move to the extended Amati relation. In Table \ref{tab:parameters_LCDM} we also show the constraints on parameters of the extended Amati relation from the low-$z$ GRB data in the low-redshift calibration method and from all GRB data in the global fitting method, which are completely consistent with each other at about $1\sigma$ confidence level. This also implies that the extended Amati relation could ease the tension of constraints on parameters of GRB relation obtained from different GRB data combinations, similar with the conclusion in the previous section \ref{sec:evolution}. Then we compare the constraints on the matter density parameter using these two methods. In Fig. \ref{fig:Omm} we plot the two-dimensional marginalized distribution with 1$\sigma$ and 2$\sigma$ contours between $\Omega_{\rm m}$ and $\alpha$, as well as the one-dimensional distributions of $\Omega_{\rm m}$ and $\alpha$ (blue lines). We also show the constraint on $\Omega_{\rm m}$ using the low-redshift calibration method (red dashed line) in the lower-right panel for comparison.

Clearly, we can see that the constraint on $\Omega_{\rm m}$ using the global fitting method is nearly unconstrained, while the powerful constraint on $\Omega_{\rm m}$ is obtained by using the low-redshift calibration method. After our careful numerical checks, we find that the reason of this huge difference is that the global fitting method includes all the correlations among parameters, such as the matter density parameter and parameters of GRB relation (as shown in the lower-left panel of Fig. \ref{fig:Omm}), during the calculations, while this low-redshift calibration method totally neglects these strong degeneracies by force and only uses the simple error propagation equation [see Eq.(\ref{error})], which only includes the standard deviations of parameters themselves, to estimate the final statistical errors for the data. This will inevitably underestimate the statistical errors of data and then overestimate the constraining power on the matter density parameter. This means that the low-redshift calibration method gives the strongly biased and overestimated constraint on the matter density parameter, which is untrustable.

In the two dimensional contour (lower left panel) of Fig. \ref{fig:Omm}, we add a vertical line which denote the best-fit value $\alpha=0.72$ obtained from the low-redshift GRB data alone in the low-redshift calibration method. This line crosses the $1\sigma$ and $2\sigma$ contours with the crossing values 0.11 $\leq\Omega_{\rm m}\leq$ 0.36 and 0.05 $\leq\Omega_{\rm m}\leq$ 0.76. We find that these values are quite similar with the $1\sigma$ and $2\sigma$, lower and upper limits of the constraint on the matter density parameter by using the low-redshift calibration method (Eq.~\ref{eq:omm:low-z}). This analysis implies that in the low-redshift calibration method, the strong correlations between $\Omega_{\rm m}$ and coefficients of the extended Amati correlation are not fully taken into account in the calculation. Neglecting this information could break the strong degeneracies between $\Omega_{\rm m}$ and correlations coefficients by force and then give quite stringent constraint on the matter density parameter.

Furthermore, in the low-redshift calibration method, we find that the final constraint on $\Omega_{\rm m}$ is very sensitive on how to calibrate the relation. Here, we assume three sets of parameters to calibrate the extended Amati relation:
\begin{itemize}
  \item (a) the best-fit values obtained from the low-$z$ GRB data;
  \item (b) the best-fit values obtained from all GRB data;
  \item (c) $a=52.40$, $b=1.45$, $\alpha=0.10$ and $\beta=0.14$, which is quite similar with the first case (a), except the $\alpha$.
\end{itemize}
Based on the Table \ref{tab:parameters_LCDM}, we can easily know these three cases are consistent with each other at $1\sigma$ confidence level. Now we use the low-redshift calibration method to constrain $\Omega_{\rm m}$, which is shown in the lower-right panel of Fig. \ref{fig:Omm} [(a) the red dashed line, (b) green dash dot line and (c) black dotted line, respectively]. The case (a) is the standard case in the low-redshift calibration method and obtains the stringent constraint on $\Omega_{\rm m}$ (see Eq.~\ref{eq:omm:low-z}). But if we use case (b) to calibrate the relation, the obtained constraint on the matter density parameter is quite different from that in case (a):
\begin{equation}
  \Omega_{\rm m} = 0.58\pm0.12~(68\%~{\rm C.L.})~.
\end{equation}
The small shift of the parameters of ``calibrated'' relation could significantly change the final constraint on $\Omega_{\rm m}$ in the low-redshift calibration method, which implies that the constraint on $\Omega_{\rm m}$ is very sensitive to the values of parameters for the calibration, when we do not fully take the strong correlations among parameters into account and underestimate the statistical errors. When we use the case (c), the situation becomes worse. The obtained constraint on $\Omega_{\rm m}$ is weaker and only has the lower limit: $\Omega_{\rm m} > 0.71$ at 95\% confidence level.

\begin{figure*}
	\centering
	\includegraphics[width=0.9\textwidth]{sim_distribution_z_Ep.eps}
	\caption{The distributions of redshift $z$ (left panel) and peak energy $E_{\rm p}$ (right panel). The blue histogram represent the distributions of redshift $z$ and peak energy $E_{\rm p}$ respectively. the red lines are the distribution function of this two quantities.}\label{fig:distributions}
\end{figure*}

We also perform the similar check in the original Amati relation. The constraints on the matter density parameter using the cases (a) and (b) to calibrate the relation are shown in the upper-right panel of Fig.\ref{fig:Omm} [(a) the red dashed line and (b) green dash dot line, respectively]. We obtain the same conclusion that the constraint on $\Omega_{\rm m}$ is very sensitive to the values of parameters for the calibration in the low-redshift calibration method.

\subsection{Simulation Results}\label{sec:section3.3}

For the time being, the number of GRB data we can use for possible cosmological studies are very small, when comparing with various type Ia supernovae data sets. There are still very large intrinsic scatter and large statistical errors on data in the current GRB data set, as we discuss above. In order to check how large biased constraint on the matter density parameter by using the low-redshift calibration method, we simulate the mock GRB data with high precisions.

To get the simulation data of GRBs, we first plot the distribution of the redshift $z$ and the observed quantity peak energy $E_{\rm p}$ of the 138 GRBs sample \citep{Liu:2015} in Fig. \ref{fig:distributions} with the blue histogram. Then we find two distribution functions to trace the histogram of redshift $z$ and $E_{\rm p}$. For the redshift $z$, we assumes the distribution function has the form
\begin{equation}\label{equ:distribution_z}
p(z) \propto z e^{-z}.
\end{equation}
We also assumes the distribution of peak energy $E_{\rm p}$ is gamma distribution
\begin{equation}\label{equ:gammadistribution}
p(x;\alpha,\lambda)=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\lambda x},
\end{equation}
where $\alpha$ and $\lambda$ are parameters, and the gamma function is
\begin{equation}
\Gamma (\alpha) = \int_{0}^{\infty} e^{-t}t^{\alpha-1}dt.
\end{equation}
We put the distribution function of $E_{\rm p}$ in Fig. \ref{fig:distributions} (right panel) the form of
\begin{equation}\label{equ:distribution_Ep}
p(E_{\rm p}/100) \propto \frac{1}{\Gamma(\alpha)}(E_{\rm p}/100)^{\lambda}e^{-E_{\rm p}/100}.
\end{equation}
We also use the distribution of $E_{\rm p,i}$ to generate the mock GRB data, and find that the obtained constraints on the parameters of relation and the matter density parameter are consistent with that by using the distribution of $E_{\rm p}$.

We use the fiducial values of the extended Amati relation: $a=52.52$, $b= 1.59$, $\alpha= 0.33$, $\beta= -0.31$, obtained from all-redshift GRB data by using the global fitting calibration method in the flat $\Lambda$CDM framework. (Different choices of fiducial values do not affect the final results.) Based on these two distribution functions, equations (\ref{equ:distribution_z}) and (\ref{equ:distribution_Ep}), we could obtain the simulation data for both redshift $z$ and peak energy $E_{\rm p}$, and then get the simulation data of bolometric fluence $S_{\rm bolo}$ according to Eq.~(\ref{equ:newrelation}). For the real GRB data, we find that the uncertainties of $E_{\rm p}$ and $S_{\rm bolo}$ have the mean relative error of $28\%$ and $18\%$, respectively \citep{Liu:2015}. Therefore, here we set the simulation data of $E_{\rm p}$ and $S_{\rm bolo}$ have the same relative errors. Furthermore, in the process of simulation, we also include the information of the intrinsic scatter and choose $\sigma_{\rm int}=0.353$ obtained from the calibration method using all current GRB data. Finally, the $\chi^2$ equation becomes
\begin{equation}\label{equ:sim_chi2}
\chi^{2}(a,b,\alpha,\beta)=\sum_{i}\frac{\left(y_{i}-A-Bx_{i}\right)^{2}}{\sigma^{\prime 2}_{y_{i}}+B^{2}\sigma_{x_{i}}^{2}}.
\end{equation}
where $\sigma^{\prime}_{y_{i}}=\sqrt{\sigma^{2}_{\rm int}+\sigma^{2}_{y_{i}}}$ which includes the intrinsic scatter.

Firstly, we use the current precision ($\sigma_{E_{\rm p}}$,$\sigma_{S_{\rm bolo}}^{\prime}$) to simulate the 138 GRB data to constrain the coefficients of the extended Amati relation and obtain the similar best-fit values and error bars of coefficients with those from the current real GRB data, which means that we could safely use this mock data for further calculations.

Next, we use the low-redshift calibration method to constrain the matter density parameter from this mock GRB data. In practice, we select the low-redshift GRB samples with $z<1.4$ to constrain those four coefficients of the extended Amati relation. The results are shown in Fig. \ref{fig:1sigma_point_1}. The obtained best-fit values perfectly recover the fiducial ones of parameters of GRB relation. In section \ref{sec:section3.2}, we have proven that a small shift of the parameters of ``calibrated'' relation could significantly change the final constraint on $\Omega_{\rm m}$ from the real GRB data in the low-redshift calibration method. Therefore, here we want to check this using the mock GRB data. We assume the values of parameters (52.805, 2.408, -0.020, -1.939), which is a small shift from the best fit fiducial values and safely lay in the $1\sigma$ contour (red points in Fig. \ref{fig:1sigma_point_1}), to calibrate the extended Amati relation. Then we use the mock high-redshift GRB data to constrain the matter density parameter, based on this ``calibrated'' relation, and obtain the 68\% C.L. constraint:
\begin{equation}
\Omega_{\rm m}=0.65\pm0.15~,
\end{equation}
which is far away from the fiducial value $\Omega_{\rm m}=0.308$. This result tells us that even a small shift of coefficients in the $1\sigma$ contour, the constraint on the matter density parameter with the current precision of GRB data could still be significantly biased using the low-redshift calibration method. The constraint on $\Omega_{\rm m}$ is very sensitive to the values of parameters for the calibration, when we do not fully take the strong correlations among parameters into account and underestimate the statistical errors in the low-redshift calibration method.

In Fig \ref{fig:find_min_grb_all} we show the obtained biased values of matter density parameter from the mock GRB data with different precisions using the low-redshift calibration method. We can see that when GRB data have larger samples with the current precision or the relative errors of GRB data decreases, $\Omega_{\rm m}$ will be more and more close to the fiducial value.
\begin{figure}
	\centering
	\includegraphics[width=0.45\textwidth]{1sigma_point_138.eps}\\
	\caption{The 1-D and 2-D marginalized distributions with $1\sigma$ and $2\sigma$ contours for the coefficients $a$, $b$, $\alpha$ and $\beta$ from the mock GRB data with the current precision of GRB data.}\label{fig:1sigma_point_1}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=0.45\textwidth]{find_min_grb_all.eps}\\
	\caption{The best-fit value and 1$\sigma$ error of $\Omega_{\rm m}$ for different quantities of GRBs. The red points with error bars represent the best-fit value and 1$\sigma$ error of $\Omega_{\rm m}$ for different numbers of GRBs respectively. The red dashed line corresponding to the fiducial value of $\Omega_{\rm m}=0.308$.}\label{fig:find_min_grb_all}
\end{figure}

\section{Conclusions and discussions}\label{sec:summary}

Thanks to the extremely high power of the explosion, Gamma-ray bursts (GRBs) are proposed as a promising candidate to trace the Hubble diagram of the Universe in high redshift range. In recent years, several popular luminosity correlations have been statistically concluded from GRBs observations. People also have made great efforts to standardize them as cosmological distance indicators. In general, there are two key issues extensively discussed in the literature when we use these correlations for cosmology. The first issue is that whether luminosity correlations evolve with redshift. The other one is about methods to calibrate luminosity correlations. In this paper we take the Amati relation as an example and use the current GRB data in the redshift range [0.0331,~8.1] to investigate these two issues. Here we summarize our main conclusions in more detail:
\begin{itemize}
  \item We divide the whole GRB data into two redshift bins and constrain the coefficients $a$ and $b$ of the original Amati relation in each bin, respectively. Based on the MCMC method, we find that the constraints on both the intercept $a$ and the slope $b$ are quite different from the low-redshift and high-redshift GRB data, respectively. The tension is more than $1\sigma$ level, which implies that the parameters of original Amati relation could be redshift-dependent.
  \item We introduce two extra redshift-dependent terms to characterize the redshift evolution of the luminosity correlation. Interestingly, we find that the tension between constraints of parameters from low-$z$ and high-$z$ GRB data are alleviated.
  \item Besides the evolution of Amati correlation, the calibration is also very important for using GRB data in cosmological studies. We firstly check the constraint on the matter density parameter using the low-redshift calibration method and obtain $\Omega_{\rm m}=0.24\pm0.12$ (68\% C.L.) which is consistent with other works.
  \item However, we could also use the global fitting method to calibrate the GRB data, in which we use all the GRB data to constrain coefficients of the extended Amati relation and the cosmological parameters simultaneously. In this case, we find that, due to the current poor precision of the GRB data, the constraint on $\Omega_{\rm m}$ is very weak, which is quite different from that by using the low-redshift calibration method.
  \item After our careful checks, we find that the low-redshift calibration method does not take the whole correlations between $\Omega_{\rm m}$ and coefficients into account. Neglecting the correlation information can break the degeneracies between $\Omega_{\rm m}$ and coefficients. Therefore, the obtained constraint on $\Omega_{\rm m}$ is totally biased. A small shift of the parameters of ``calibrated'' relation could significantly change the final constraint on $\Omega_{\rm m}$ in the low-redshift calibration method, which implies that the constraint on $\Omega_{\rm m}$ is very sensitive to the values of parameters for the calibration.
  \item In order to investigate how large biased constraint on $\Omega_{\rm m}$ by using the low-redshift calibration, we simulate the mock GRB data with different precisions. We find that the mock data with current precision will give significantly biased result using the low-redshift calibration method, even we assume a small shift on the parameters of ``calibrated'' relation. When GRB data include larger samples with the current precision or the relative errors of GRB data decreases, the constraint on $\Omega_{\rm m}$ is more and more close to the fiducial value.
\end{itemize}

\section*{Acknowledgements}
G.-J. Wang and Z.-H. Zhu are supported by the the National Science Foundation of China under grant No. 11373014. H. Yu is supported by the National Basic Research Program of China (973 Program, grant No. 2014CB845800). Z.-X. Li is supported by the NSFC under grant No. 11505008. J.-Q. Xia is supported by the National Youth Thousand Talents Program and the NSFC under grant No. 11422323. The research is also supported by the Strategic Priority Research Program ``The Emergence of Cosmological Structures'' of the Chinese Academy of Sciences, Grant No.XDB09000000, and the NSFC under grant Nos. 11633001 and 11690023.


\section{Introduction}
The anomalous Hall effect was first discovered by Hall in 1881, and occurs as a result of broken time-reversal symmetry in metallic
ferromagnets with strong spin-orbit coupling\cite{Nagaosa2010}. The time-reversal
symmetry can be broken by the inclusion of an appropriate weak magnetic field\cite{Zhai2005}.
In the early stages of research on this phenomenon, "anomalous behavior" was reported by Karplus and Luttinger\cite{Karplus1954}, which naturally arose in the first microscopic theory of the anomalous Hall effect.
The Berry curvature $\Omega_{n}(\mathbf{k})$  concept was adopted theoretically to relate the anomalous Hall effect to the topological nature of the Hall currents. As a result of the topological Berry phase, the intrinsic anomalous Hall effect can be expressed as an integral over the Fermi surface\cite{Haldane2004}.
Further, in ferromagnets, electrons have an "anomalous velocity" perpendicular to the electric field due to their $\Omega_{n}(\mathbf{k})$. In previous studies, Culcer et al.\cite{Culcer2003} have investigated the anomalous Hall effect in paramagnetic two-dimensional
systems with a external magnetic field, while Wang et al.\cite{Wang2014} have investigated the quantum anomalous
Hall effect in conventional diluted magnetic semiconductors with magnetically doped InAs/GaSb
quantum wells based on the Kane model.
Relevant background information and the basic concept of $\Omega_{n}(\mathbf{k})$ are presented in the review compiled by Xiao et al. \cite{Xiao2010}.

All of the studies mentioned above were performed magnetically. However, the Berry phase can
manifest in various magnetic and nonmagnetic materials. As regards cases without a magnetic field, the Berry phase manifests in gapped Dirac materials, where the time reversal asymmetry is realized via a nonequilibrium
valley polarization\cite{Mak2014,Lensky2015}. Thus, gapped Dirac materials with several valleys, such as
graphene and transition metal dichalcogenide monolayers, have been investigated by Xiao et al.\cite{Xiao2007,Xiao2012}. The authors found a valley dependent Berry phase effect which can
yield a valley contrasting Hall transport by an electric mean\cite{Xiao2007}.
Further, Kumar et al.\cite{Kumar2016} have observed chirality in bulk and edge plasmons
without an external magnetic field, where the valley imbalance induces a net $\Omega_{n}(\mathbf{k})$ (Berry flux $F$).
Note that $F$ can cause chiral plasmons
in the absence of a magnetic field\cite{Kumar2016,Song2016}. Finally, Song et al. have confirmed the existence of counterpropagating charge density waves in chiral plasmons, where splitting dispersion occurs
for transverse edge modes in opposite directions. These modes are regarded as chiral Berry plasmons (CBPs)\cite{Song2016}.

However, quantum effects may yield other interesting phenomena in these CBPs.
In this Letter, we employ a quantum hydrodynamic (QHD) method coupled with the anomalous velocity to examine the quantum effect on the splitting dispersive relation and the transverse mode confinement of CBPs propagating along the boundary of a two-dimensional electron gas (2DEG) confined to a half-plane. Note that Haas et al. first introduced the QHD theory by solving the nonlinear Schr\"{o}rdinger-Poisson or Wigner-Poisson kinetic models\cite{Haas2000,Manfredi2001}. Further, a quantum effect always exists in a metallic electron gas\cite{Haas2011}, establishing new chiral edge modes and enhancing both the splitting and the transverse confinement significantly for a wide range of transverse modes.
Finally, Gauss units will be adopted throughout this paper,
except in the case of specific definitions.

We consider CBPs along the edges of 2DEG while examining the quantum effect without a magnetic field.
The excitations of 2DEG
are described by the QHD model\cite{Haas2000,Manfredi2001} using
the continuity equation
\begin{equation}
\frac{\partial n_{e}}{\partial t}+\nabla\cdot
(n_{e}\mathbf{V})=0,  \label{ne
\end{equation}
the momentum-balance equation
\begin{equation}
\frac{\partial \mathbf{u}_{e}}{\partial t}+(\mathbf{u}_{e}\cdot
\nabla) \mathbf{u}_{e}=\frac{e}{m_{e}}\nabla\phi
-\frac{\pi \hbar ^{2}}{m_{e}^{2}}\nabla n_{e}
-\frac{\hbar^{2}}{2m_{e}}\nabla (\frac{1}{\sqrt{n_{e}}}\nabla^{2}\sqrt{n_{e}})
,  \label{mom}
\end{equation
and the integral form of the electric potential
\begin{equation}
\phi(\mathbf{r},t)=\int d^{2}\mathbf{r}' W(\mathbf{r}-\mathbf{r}')(n_{e}(\mathbf{r'},t)-n_{0}).
\end{equation}
Here, $m$ is the effective electron mass, $e$ is the elementary charge, $\hbar$ is
the Planck constant, $W(\mathbf{r}-\mathbf{r}')$ is the Coulomb interaction,
$n_{e}$ is the density, and $\mathbf{u}_{e}$ is the homogeneous fluid velocity.
In particular,
\begin{equation}
\mathbf{V}=\mathbf{u}_{e}+\mathbf{V}_{a}, \mathbf{V}_{a}=\frac{eF}{\hbar}[(\nabla\phi)\times \mathbf{e}_{z}],
\end{equation}
with chiral velocity $\mathbf{V}_{a}$, which is a self-induced anomalous
velocity component\cite{Xiao2010}, and $F$ is dimensionless\cite{Song2016}.
Note that the final two terms on the right hand side of Eq. (\ref{mom}) are regarded as quantum effects\cite{Haas2000,Manfredi2001}. That is, the second term is the quantum statistical effect (second-order quantum correction), which is the force due to the internal interactions in the electron species, and the third term is the quantum diffraction effect (fourth-order quantum correction) due to the quantum pressure (the so- called "Bohm potential").
For definiteness, only the quantum statistical effect (second order) is considered in this semi-analytical study, as the inclusion of the quantum pressure term (fourth order) cannot yield an analytical solution.
The above nonlinear equations can be linearized under weak perturbation of the plasmons, where $n_{e}(\mathbf{r},t)$ is represented by the first-order
perturbed value $n_{e}=n_{0}+n_{e1}$, with $n_{e1}\ll n_{0}$\cite{Fetter1985}. $n_{0}$ is the initial equilibrium density.
Thus, the linearized QHD theory yields
\begin{equation}
\frac{\partial n_{e1}}{\partial t}+n_{0}\nabla\cdot
\mathbf{V}=0,  \label{num2d
\end{equation}
\begin{equation}
\frac{\partial \mathbf{u}_{e}}{\partial t}=\frac{e}{m}\mathbf{\nabla}\phi
-\frac{\pi \hbar ^{2}}{m^{2}}\nabla n_{e1}
,  \label{mom2d}
\end{equation
and
\begin{equation}
\phi(\mathbf{r},t)=\int d^{2}\mathbf{r}' W(\mathbf{r}-\mathbf{r}')n_{e1}(\mathbf{r'},t). \label{phi}
\end{equation}

A semi-infinite ($x\geq0$) 2DEG in the half $x$-$y$ plane is considered, where
 $n_{e}(\mathbf{r},t)$ and $\mathbf{V}(\mathbf{r},t)$ take finite values for  $x\geq0$, while $n_{e}(\mathbf{r},t)=0$ and $\mathbf{V}(\mathbf{r},t)=0$ for $x<0$, as shown in Fig. 1. Thus, an edge plasmon is formed along the boundary at $x=0$, propagating as a plane wave along the $y$ axis with frequency $\omega$ and wave number $q$. The potential, perturbed density, and velocities can be expressed as
\begin{equation}
\begin{array}{lcl}
    \phi(\mathbf{r},t)=\phi_{q}(x)\exp(i\omega t-iqy), \\
    n_{e1}(\mathbf{r},t)=n_{e1q}(x)\exp(i\omega t-iqy), \\
    u_{e}(\mathbf{r},t)=u_{eq}(x)\exp(i\omega t-iqy), \\
    \mathbf{V}(\mathbf{r},t)=\mathbf{V}_{q}(x)\exp(i\omega t-iqy),
 \end{array}
 \label{Fourier}

\end{equation
respectively. Thus, Eq. (\ref{num2d}) can be rewritten as $n_{e1}(\mathbf{r})=-n_{0}\mathbf{\nabla}\cdot\mathbf{V}(\mathbf{r},t)/i\omega$. Further, using the step-like property of $\mathbf{V}(\mathbf{r},t)$ [$\mathbf{V}(\mathbf{r},t)\neq0$ for $x\geq0$ and $\mathbf{V}(\mathbf{r},t)=0$ for $x<0$], a jump condition for the $x$ derivative of $\phi_{q}(x)$ can be obtained based on  Eq. (\ref{phi})\cite{Song2016}, where
\begin{equation}
\frac{\partial \phi_{q}}{\partial x}\mid_{x=0^{+}}-\frac{\partial \phi_{q}}{\partial x}\mid_{x=0^{-}}=\frac{1}{i\omega}(\frac{\partial W_{q}}{\partial x}\mid_{x=0^{-}}-\frac{\partial W_{q}}{\partial x}\mid_{x=0^{+}})n_{0}\mathbf{V}_{x}\mid_{x=0^{+}}. \label{refe7
\end{equation}
Here, $\mathbf{V}_{x}\mid_{x=0^{+}}=\mathbf{u}_{ex}\mid_{x=0^{+}}+\mathbf{V}_{ax}\mid_{x=0^{+}}$ and $W_{q}=-e\int dk exp(ikx)/\sqrt{q^{2}+k^{2}}$, assuming a wave number $k$ along the $x$ axis.

According to the literature\cite{Song2016,Fetter1985}, $W_{q}$ can be simplified as
\begin{equation}
W_{q}(x)\approx-4\pi e\int\frac{dk}{2\pi}\frac{|q|}{k^{2}+2q^{2}}\exp(ikx)=-\frac{4\pi e}{2\sqrt{2}}\exp(-|q|x), x\geq0,
\end{equation}
for small $ k/q $.
Again, by substituting the simplified $W_{q}(x)$ into the integral potential, we find
\begin{equation}
\phi_{q}(x)=\int dx' W_{q}(x-x')n_{e1q}(x').
\end{equation}
Then, using the convolution theorem, we can obtain a differential
equation of $\phi_{q}(x)$ satisfying
\begin{equation}
\begin{array}{lcl}
 (x^{2}-2q^{2})\phi_{q}(x)=0, x<0, \\
 (x^{2}-2q^{2})\phi_{q}(x)=4\pi e |q|n_{e1q}(x), x\geq0.
 \end{array}
 \label{phi12}

\end{equation
Thus, $\phi_{q}(x)$ has solutions of
\begin{equation}
\begin{array}{lcl}
 \phi_{q}(x)=\phi_{1}\exp(\kappa_{1}x), x<0,\\
\phi_{q}(x)=\phi_{2}\exp(-\kappa_{2}x), x\geq0.
 \end{array}
 \label{phi1212}

\end{equation
Here, $\kappa_{1}=\sqrt{2}|q|$
and $\kappa_{2}$ is determined as follows.
Applying $\partial/\partial t$ to the continuity relation in Eq. (\ref{num2d}) and
substituting $\frac{\partial \mathbf{V}}{\partial t}=\frac{\partial \mathbf{u_{e}}}{\partial t}+\frac{\partial \mathbf{V_{a}}}{\partial t}$ into Eq. (\ref{num2d}), we find
\begin{equation}
  \frac{\partial^{2}n_{e1}}{\partial t^{2}}=-\frac{e}{m}n_{0}\nabla
^{2}\phi+\frac{\pi \hbar^{2} n_{0}}{m^{2}}\nabla^{2}n_{e1}
. \label{refe4
\end{equation}
Using the plane-wave forms of $n_{e1}(\mathbf{r},t)$ and $\phi(\mathbf{r},t)$ and the relation $ n_{e1q}(x)=(x^{2}-2q^{2})\phi_{q}(x)/(4\pi e |q|)$ based on Eq. (\ref{phi12}), $n_{e1}(\mathbf{r},t)$ can be expressed as
\begin{equation}
n_{e1}(\mathbf{r},t)=\frac{\phi_{2}}{4\pi e|q|}(\kappa_{2}^{2}-2q^{2})e^{-\kappa_{2}x}e^{i\omega t-iqy}, \label{ne1q}
\end{equation}
for $x\geq0$.
Therefore, Eq. (\ref{refe4}) can be rewritten as
\begin{equation}
1=\frac{2\omega_{P}^{2}}{\omega^{2}}\frac{\kappa_{2}^{2}-q^{2}}{\kappa_{2}^{2}-2q^{2}}
-A_{Q}\frac{\kappa_{2}^{2}-q^{2}}{\omega^{2}}, \label{kappa2}
\end{equation}
where the bulk plasmon frequency $\omega_{P}=\sqrt{2\pi n_{0}e^{2}|q|/m}$, according to the literature\cite{Fetter1985}. Here, we assume that $q\ll k_{F}=\sqrt{2\pi n_{0}}$, where $k_{F}$ is the Fermi wave number. In particular, $A_{Q}=\pi\hbar^{2}n_{0}/m^{2}$ is the coefficient of the quantum term; therefore, $A_{Q}(\kappa_{2}^{2}-2q^{2})$ characterizes the quantum effect term.
Equation (\ref{kappa2}) yields a nonlinear equation of $\kappa_{2}$. Further, $\kappa_{2}^{2}$ can be explicitly expressed as
\begin{equation}
\kappa_{2}^{2}=\frac{3A_{Q}q^2 - \omega^2 + 2 \omega_{P}^2 \mp
  \sqrt{A_{Q}^2q^4 + 2A_{Q}q^2\omega^2 + \omega^4 + 4A_{Q}q^2\omega_{P}^2 - 4\omega^2\omega_{P}^2 + 4\omega_{P}^4}}{2A_{Q}}. \label{kappa22}
\end{equation}
 At $q=0$, $\kappa_{2}^{2}$ has two solutions: $\kappa_{2}^{2}=0$ and $\kappa_{2}^{2}=-2\omega^{2}$. We consider positive $\kappa_{2}^{2}$ only, as the mode for negative $\kappa_{2}^{2}$ joins the bulk. Thus, the potential decays exponentially from the boundary $x=0$, yielding non-vanishing edge modes.

In what follows, in order to obtain the dispersion relation, we first express
\begin{equation}
{V}_{x}|_{x=0^{+}}={u}_{ex}|_{x=0^{+}}+{V}_{ax}|_{x=0^{+}}.
\end{equation}
Here,
\begin{align}
{u}_{ex}|_{x=0^{+}}=\frac{1}{i\omega}[-\frac{e}{m}\kappa_{2}\phi_{2}
+\frac{\pi\hbar^{2}}{m^{2}}\frac{1}{4\pi e|q|}\kappa_{2}(\kappa_{2}^{2}-2q^{2})\phi_{2}], \label{ue}
\end{align}
based on the Fourier transform of the momentum-balance equation and Eq. (\ref{ne1q}),
and
\begin{equation}
{V}_{ax}|_{x=0^{+}}=\frac{-iqeF}{\hbar}\phi_{2}.
\end{equation}
The dispersion relation is determined from the continuity condition of $\phi$ and the jump boundary condition given in Eq. (\ref{refe7}), such that
\begin{equation}
\sqrt{2}|q|+\kappa_{2}=\frac{2\omega_{p}^{2}}{\omega^{2}}\kappa_{2}-A_{F}\frac{q^{2}sgn(q)}{ \omega}
-A_{Q}\frac{\kappa_{2}(\kappa_{2}^{2}-2q^{2})}{\omega^{2}}, \label{wq}
\end{equation}
indicating a nonlinear equation of $\kappa_{2}$ and $\omega$. Here, $A_{F}=4\pi e^2F/\hbar$ is the coefficient of the chiral Berry term.
At $q=0$, the above dispersion relation yields
$\kappa_{2}=-A_{Q}\kappa_{2}^{3}/\omega^2$, indicating a diverging mode or $\kappa_{2}=0$, which will induce a jump in the edge-mode frequency at $q=0$.
We can couple the two nonlinear equations (Eqs. (\ref{kappa2}) and (\ref{wq})) in order to semi-analytically solve $\kappa_{2}$ and $\omega$ for a given $q$. The solutions obey two coupled polynomial functions of the fifth and fourth orders, respectively.
It is instructive to discuss the properties of Eqs.(\ref{kappa2}) and (\ref{wq}) in the vicinity of $q=0$.
Indeed, the coupling of these two equations yields the solution
\begin{align}
\omega=
&- (2.8A_{Q}q\kappa_{2}^4)/A_{F} - (2.8A_{Q}q\kappa_{2}^3)/A_{F} + (0.5(- 1.4A_{F}^2q^2 + 9.9A_{Q}q^2 + 5.7\omega_{P}^{2})\kappa_{2}^2)/(A_{F}q) \notag \\
&+ (2.9(A_{Q}q^2 + \omega_{P}^{2})\kappa_{2})/(A_{F}q) - (- 0.4A_{F}^2q^2 + 2.1A_{Q}q^2 + 2.1\omega_{P}^{2})/(A_{F}q). \label{wqe}
\end{align}
Here, $\omega$ represents $\omega_{+edge}$ and $\omega_{-edge}$ for negative and positive $q$, respectively.

Thus, we can obtain the limits of both $\omega$ and $\kappa_{2}$ for $q\rightarrow 0^{-}$ and $q\rightarrow 0^{+}$:
\begin{equation}
\begin{array}{lcl}
 \kappa_{2}|q\rightarrow 0^{-}=O(q^{4}),   \\
\omega|q\rightarrow 0^{-}=\frac{21\omega_{p}^2}{10A_{F}|q|},  \\
 \kappa_{2}|q\rightarrow 0^{+}=\frac{\sqrt{2}}{2}q,  \\
\omega|q\rightarrow 0^{+}=\frac{\sqrt{2}}{8}A_{F}q.
 \end{array}
 \label{limits}

\end{equation
The results presented below are determined by taking
$n_{0}=6\times10^{10}$ cm$^{-2}$ and $m=0.03m_{e}$, based on the literature\cite{Yan2012,Song2016}, where $m_{e}$ is the free electron mass.
For definiteness, we take the unit of $q$ as $k_{F}$ for large $q$, but as cm$^{-1}$ for small $q$, and the unit of the plasmon frequency $\omega$ is $E_{F}/\hbar$ with the Fermi energy $E_{F}=\hbar^{2}k_{F}^2/2m$. Thus, the strength of the quantum effect is defined by a dimensionless parameter $\widetilde{Q}=A_{Q}(\hbar k_{F}/E_{F})^2$.

Figures 2(a, b) show the dispersion relations of CBPs with $F=1$, comparing the edge modes obtained without (circles) and with (triangles) the quantum effect.
The $\omega_{+ edge}$ mode without the quantum effect diverges between $- 0.0028<q/k_{F}<-0.0013$, and both $\omega_{\pm edge}$ modes with the quantum effect diverge between $0.0024<|q|/k_{F}\leq0.0033$ and $q=0$.
Indeed, with the quantum effect, substituting these $q$ values into
Eqs. (\ref{kappa2}) and (\ref{wq}) yields negative decay constant $\kappa_{2}$ and edge frequencies $\omega_{\pm edge}$, inducing the diverging modes.
Without the quantum effect, in order to validate our model, we first compare our results with those of the literature\cite{Song2016} for small $q$. The splitting trends of the edge modes agree well with the linear approximation results (the first-order $q$ modes) given in Ref.\cite{Song2016}, as shown in Fig. 2(a) (circles and squares ($\omega_{P}$)).
For large amplitude and negative $q$, a new edge mode emerges for $|q|/k_{F}>0.0028$ (circles) (Fig. 2(a)), due to the dominance of the quadratic ($\propto q^2$) terms; this behavior did not occur for the linear approximation limit.
With the quantum effect,
$\omega_{+edge}$ exhibits a maximum value in the vicinity of $q=0^{-}$,  and decreases with increasing $|q|$ for $|q|/k_{F}\leq0.0024$ (Fig. 2(a)). This behavior occurs because the dispersion relation with the quantum effect has a vanishing solution at $q=0$ (see the description below Eq. (\ref{wq})).
This scenario is reversed for the $\omega_{-edge}$ mode, where $\omega_{-edge}$ exhibits minimum values near $q=0^{+}$ and increases quickly with increasing $|q|$, for large $q$.

Figure 2(b) indicates that $\omega_{\pm edge}$, respectively, exhibits a new one-way edge state with the quantum effect compared to that without the quantum effect, at the both side of $q=0$.
Although this one-way edge CBP state discussed here presents some similar characters with
the one-way edge magnetoplasmon state reported in Ref.\cite{Jin2016}, there are fundamental differences between them. The one-way edge magnetoplasmon has been realized by breaking the time-reversal symmetry with oppositely-directed magnetic fields, whereas the one-way edge CBP is realized by considering the quantum effect.
Indeed, the $\omega_{\pm edge}$ with the quantum effect have different limits (see Eq. (\ref{limits})) for $q\rightarrow 0^{-}$ and $q\rightarrow 0^{+}$, whereas the $\omega_{\pm edge}$ without the quantum effect have symmetric distribution at the different sides of $q=0$ and tend to zero at $q=0$ (see (Inset) in 2(b)).
Note that the splittings have crucial applications in subwavelength optical
nonreciprocity\cite{Sounas2013} in the terahertz (THz) range.

The corresponding transverse confinements of these edge modes are displayed in Fig. 3 , for (a) negative and (b) positive $q$, for the same conditions as in Fig. 2. Without the quantum effect, Fig. 3(a) again shows divergence of $\omega_{+edge}$  for $- 0.0028<q/k_{F}<-0.0013$, corresponding to Fig. 2(a) (circles).
With the quantum effect, both $\omega_{\pm edge}$ modes exhibit greater localization over a wider range of $q$; this behavior is ascribed to the coupling of the quantum effect and the chiral velocity (see Eq. (\ref{wqe})). This significantly enhanced transverse localization is preferable for the realization of subwavelength optical non-reciprocal devices in the THz range with no magnetic field. Note that this possibility avoids the current requirement for large magnetic-based devices to realize nonreciprocity in the THz range, enhancing industrial potential as regards microwave and photonic components.

Without the quantum effect, the splitting between $\omega_{\pm edge}$  increases linearly with increasing $q$ and $F$, because of the relations $\omega \propto Fq$ and $\mathbf{V}_{a}=\frac{eF}{\hbar}[(\nabla\phi)\times \mathbf{e}_{z}]$, which are given in Ref.\cite{Song2016}. However, with the quantum effect, the splitting  between $\omega_{\pm edge}$ indicates a nonmonotonous relation to $q$ and $F$, because of the diverging solutions at certain $q$  and $F$ values. Indeed, diverging solutions usually occur for a higher-order polynomial function. Fig. 4  shows comparisons of (a) the edge modes $\omega_{\pm edge}$ and (b) the corresponding transverse confinement lengths, with and without the quantum effect, as a function of $F$, and for fixed $q/k_{F}=0.0098$.
Without the quantum effect,
in Fig. 4(a), the splitting modes (squares and spheres) agree well with those in the literature\cite{Song2016}. However, a new $\omega_{+ edge}$ mode branch emerges as a result of the quadratic $q^{2}$ terms for $F\geq 0.55$; this new mode was easily omitted from the linear $q$ approximation. As a result, this new mode introduces a new confinement length branch for $F\geq 0.55$, as can be clearly seen in Fig. 4(b) (squares). In contrast, the confinement lengths of the $\omega_{- edge}$ mode and the increasing branch of the $\omega_{+ edge}$ are in good agreement with the curves in Ref.\cite{Song2016}. This increasing branch diverges when $F\geq 0.35$, and a new $\omega_{+ edge}$ mode branch develops when $F\geq 0.55$. The latter is more closely compressed to the edge than its counterpropagating $\omega_{- edge}$ mode.
With the quantum effect, $\omega_{+ edge}$ vanishes for $F< 0.55$. This behavior occurs because small $F$ yields negative decay constants and vanishing frequencies. Again, the oscillations in the $\omega_{- edge}$ mode can be attributed to the fact that the dispersion relation is determined by polynomial functions of the high-order decay constant $\kappa_{2}$, where the coefficients change non-monotonically with increasing $F$,  as shown in Eq.(\ref{wqe}). With increasing $F$, the $\omega_{+ edge}$ mode introduces a cut-off confinement of 0.4 (upright triangles), and this mode is primarily confined to the edge in both cases, i.e., both with and without the quantum effect. Such a confinement will have crucial importance for the realization of nonreciprocal devices with no magnetic field.

Finally, we briefly discuss the properties of CBPs with the quantum effect, compared to magnetoplasmons.
For magnetoplasmons, the authors of Ref.\cite{Fetter1985} have reported that the quantum effect with strength defined as being $\propto sq$ has almost no effect on the edge modes in the limit $sq<<\omega_{P}$, with $s^2=m^{-1}(\partial p/\partial n)$. This means that the edge magnetographs are almost independent of the quantum effect.
In contrast, for the CBPs, the quantum effect induces several new edge modes (see the $\omega_{\pm edge}$ obtained with the quantum effect in the above figures), even for the limit $sq<<\omega_{P}$ (defined as $A_{Q}q^2<<\omega_{P}^2$ in the present work). Indeed, the quantum effect clearly enhances the chirality as well as the transverse confinement, indicating strong coupling of both these properties (see Eq. (\ref{wqe})).
The chirality in CBPs emerges from the coupling action of the self-induced electric field $-\nabla \phi$, due to the collective density oscillation and the anomalous velocity $\mathbf{V}_{a}=\frac{eF}{\hbar}[(\nabla\phi)\times \mathbf{e}_{z}]$, the strength of which depends on $-\nabla \phi$  and $F$. $\mathbf{V}_{a}$ can be realized by $\Omega_{n}(\mathbf{k})$ , which is an intrinsic quantum-mechanical property of a perfect crystal, because it depends on the Bloch wave function only.
Indeed, $\mathbf{V}_{a}$ can be expressed in terms of $\Omega_{n}(\mathbf{k})$ in a given Bloch state $u_{n}(\mathbf{k})$ in the $n^{th}$ band\cite{Xiao2010}, such that
\begin{equation}
\mathbf{V}_{a}=\frac{e}{\hbar} \nabla\phi\times\Omega_{n}(\mathbf{k}), \Omega_{n}(\mathbf{k})=i<\nabla_{\mathbf{k}} u_{n}(\mathbf{k})|\times|\nabla_{\mathbf{k}} u_{n}(\mathbf{k})>. \label{BC}
\end{equation}

In summary, we used a quantum hydrodynamic theory to study the chiral Berry plasmons in a semi-infinite 2DEG, with the chiral Berry term and the quantum term being considered simultaneously.
Note that this semi-infinite 2DEG can be realized easily in experiment using 2D gapped Dirac
materials, such as graphene and transition metal dichalcogenide monolayers\cite{Xiao2007,Xiao2012}.
Herein, the transverse decay constant and the edge plasmon frequency were obtained by nonlinearly and semi-analytically solving two coupled transcendental equations. First, our model was validated by comparing our results with those based on a linear approximation method with no quantum effect, for small $q$.
In our calculation, without the quantum effect, a new chiral edge mode emerged, which was due to the quadratic term contribution ($\propto q^2$).
With the quantum effect, our results predicted that both the chiral edge modes and the transverse confinement in a bounded 2DEG can be tuned. It was found that the quantum effect introduces new and intriguing phenomena. That is, a new one-way chiral edge state with the quantum effect compared to that without the quantum effect, at the both side of $q=0$ was found. Indeed, the $\omega_{\pm edge}$ modes exhibit different limits for $q\rightarrow 0^{-}$ and $q\rightarrow 0^{+}$, in contrast with the case without the quantum effect and $\omega_{P}$.
The newly discovered one-way chiral edge modes will elucidate new
paradigms in condensed matter and plasma-optics physics\cite{Sounas2013,Barnes2003,Yao2008} for the study of bounded 2DEGs, such as graphene and traditional quantum dots.

This work was supported by NSFC (11405067, 11375163, 11575135).