text:
\section{ Introduction}
  The Pancharatnam-Berry phase or geometric phase (GP) \cite{berry1984quantal,pancharatnam1956generalized}
in quantum systems  is
a natural emergence when the system is subjected to local periodic modulations of its structural parameters. In open quantum systems too,
 parametrization of the system's parameters such as tunneling coefficients \cite{kaestner2015non}, applied potential \cite{switkes}
as well as temperatures \cite{watanabe2014non} have
been known to affect open quantum systems' observables  \cite{pluecker2017gauge}.
Periodically driven temperatures in an anharmonic junction
have been shown to lead to fractional quantization
     of the heat flux \cite{ren2010berry}. Such non-trivial thermodynamics arising due to periodic modulation of temperatures
     have also been studied formally in classical Brownian heat engines where effects on
     entropy and Onsager's coefficients
    were investigated \cite{PhysRevX.5.031019} as well as in a nonequilibrium spin boson model\cite{PhysRevA.95.023610}.
    Recently,
    quantum mechanical versions of heat engines where the working medium is discretized \cite{scovil1959three} have
    resurfaced\cite{quan2007QHE,kosloff2013quantum,gelbwaser2017thermodynamic}.
The  interest in quantum heat engines (QHE)  has a twofold directive, firstly
 to answer the basic question,
 when the heat engine is really quantum \cite{PhysRevX.5.031044} and secondly to understand the thermodynamics
when the quantum effects like coherences and entanglement are predominantly active in the engine \cite{scully2003extracting,4qheEnta,gelbwaser2015strongly}.
Several intriguing works exist, such as extracting work from a single bath, increase in power aided by noise-induced quantum
coherence \cite{scully2003extracting,ScullyPNAS13092011}, surpassing Carnot efficiency using
squeezed versions of thermal baths
\cite{rossnagel2014nanoscale} as well as experimental realization of a Paul trapped single $Ca^+$ ion heat engine \cite{rossnagel2016single}. QHEs have also been experimentally
realized in  laser cooled Rb atoms\cite{PhysRevLett.119.050602} using principles based on electromagnetically induced transparency\cite{EITQHEharris}.
Despite the progress, the issue of GP in QHEs have not yet been addressed, both within and
    beyond the linear response regime of finite-time thermodynamics.
In this work we focus on a QHE where one can realize
Pancharatnam-Berry phase-like (PBp) contributions through periodically modulating the temperature of the thermal baths.

A heat engine's performance is analyzed by evaluating its efficiency and power.
Efficiency, in the finite power regime, needs to be calculated by maximizing the power with respect to some parameter of the system. This is known as the efficiency
 at maximum power ($\eta_*$), theoretically given by, $\eta_*=1-\sqrt{1-\eta_c}$\cite{Curzon-Ahlborn, Novikov1957}, where $\eta_c=1-T_c/T_h$ is the Carnot efficiency and
 $T_c(T_h)$ represents cold (hot) thermal bath's temperature.
Close to equilibrium in the endo-reversible regime, $\eta_*=\eta_c/2+\eta_c^2/8$, with the coefficient
$1/2$ being universal \cite{seifert2012stochastic,PhysRevLett.102.130602,BroeckEMP}.
  Further it is also known that
  the efficiency at maximum power is bounded above and below by
  $\eta_c/2\le\eta_*\le \eta_c/(2-\eta_c)$ \cite{izumida2012efficiency,PhysRevE.87.012133}.
  In this work, we show that  this strong standing  linear expansion
 coefficient, $1/2$, doesn't hold due to the emergence of PBp effects and one can
  go both above and below the upper bound on $\eta_*$ by having phase different driving protocols.
  Further, a QHE works by absorbing heat from the reservoirs in the form of quanta
which  mimics boson exchange, making them ideal to study heat transfer as a quantum transport phenomena from which the role of random fluctuations are studied.
Random fluctuations affect the transport properties and are usually studied via a full counting statistics (FCS)
\cite{Levitov,nazarov} method that involves calculation of moments and cumulants of $P(q,t)$, the probability distribution function (PDF) for the
  number ($q$) of particles exchanged
  between system and reservoir in a measurement time $t$.
FCS led to the steady state fluctuation theorems (FT) \cite{
seifert2012stochastic,campisi2011colloquium,uhrmp,SaarQHEPhysRevA.86.043843,UHeplQHE,campisi2015nonequilibrium}, the cornerstones of
quantum thermodynamics \cite{kosloff2013QTherm}.
         Recently, it was reported that during  transport across
  quantum junctions \cite{ren2010berry,hpg4}, the standard mathematical
  form of the FT is broken due to the emergence of PBp and hence attempts have been made to establish geometric
   FTs \cite{PhysRevE.96.022118} in spin-boson systems. We observe the same violation of the FT in our QHE. Further,  the long time PDF
  is evaluated by invoking the use of large deviation theory \cite{touchette2009large,varadhan1984large} and steadystate (SS) FTs are
  derived from large deviation results \cite{uhrmp,dhar2,lebowitz1999gallavotti}. We find that the large deviation theory cannot be used
   to determine the PDF in presence of PBp contributions.

  In this work we focus on the effects of PBp contributions on the
    thermodynamics of a coherent QHE beyond the linear response regime, where  the temperature of the two thermal baths are periodically modulated in time.
       The technique used in this paper is a standard
       procedure based on a quantum Markovian master equation (QMME) formalism (weak coupling between system and bath)
       combined with a generating function method popularly used in FCS.  The paper is organized as follows. In Sec. (\ref{qhe}), we introduce the model QHE,
       derive a QMME and discuss the PBp curvature. In Sec. ( \ref {s3}), we look at the thermodynamics of the QHE
       focusing on the geometric contributions to flux and efficiency at maximum power. In Sec. (\ref{s4}), we analyze the PBp effects on the SSFT and large deviation theory
       which is followed by conclusions and an appendix.

\section{Engine Specifications}
\label{qhe}
\begin{figure}[!tbp]
\centering
\includegraphics[width=7.5cm]{fig_0.pdf}
\caption{ A four level quantum heat
engine. Degenerate levels $|1\rangle$ and $|2\rangle$ are coupled
to two excited levels $|a\rangle$ and $|b\rangle$ through the  two driven thermal
baths  at temperatures $T_h$(t) and $T_c$(t). Levels
$|a\rangle$ and $|b\rangle$ are coupled to a cavity with a single mode
with a frequency $\nu_l$. Emission of photons into this mode is the work
done by the QHE. The PBp is quantified by the area traced in the parameter space of the two temperatures with nonzero phase difference, $\phi (\phi\ne n\pi, n\in$ integers).
}
\label{c7-schematic}
\end{figure}
The QHE model that we study has  also been studied in Refs. \citep{SaarQHEPhysRevA.86.043843,UHeplQHE,goswami2013thermodynamics,ScullyPNAS13092011}
 and is shown in Fig. (\ref{c7-schematic}).
The model consists of two thermal baths at temperatures $T_h(t)$ and $T_c(t)$. The two temperatures are periodically driven in time such that, $T_h(t)>T_c(t).$
There are two degenerate states, $|1\rangle$ and $|2\rangle$,
with forbidden transition and are coupled to the two time varying thermal baths. The higher energy state $|a\rangle (|b\rangle)$ is coupled to hot (cold) bath.
The two higher states, with an allowed transition, $|a\rangle$ and $|b\rangle$ are also coupled to a unimodal cavity.
The total Hamiltonian can be written as $\hat{H}_T= \hat{H}_o+\hat{V}$, where
\begin{equation}
\hat{H}_o=\displaystyle\sum_{\nu=1,2,a,b}
E_\nu\hat{B}_{\nu\nu}+\displaystyle\sum_{k\in h,c}\epsilon_k^{}\hat{a}_k^\dag\hat{a}
_k+\epsilon_l\hat{a}_l^\dag\hat{a}_l.
\end{equation}
Here, $\nu=1,2,a,b$, $E_\nu$ is the energy of the $\nu$th level and $\hat{B}_{\nu \nu^\prime }=|\nu \rangle\langle \nu^\prime |$ is
the  operator that causes excitation
between the states $|\nu \rangle$ and $|\nu^\prime \rangle$. $\epsilon_k$ represents the energy of $k$th mode of
the hot(h) and cold (c) thermal reservoirs and $\epsilon_l$ represents the energy of the unimodal cavity.
We model the thermal baths as harmonic modes, where
$\hat{a}^\dag (\hat{a})$ stand for the bosonic creation (annihilation) operators.
The interacting Hamiltonian is $\hat V=\hat V_{sb}^{}+\hat V_{sc}$, with the coupling between
the (working) system and thermal baths being,
\begin{equation}
\hat{V}_{sb}^{}=\displaystyle\sum_{k\in h,c}\sum_{i=1,2}\sum_{x=a,b}g_{ik}\hat{a}_k\hat{B}_{ix}^\dag+ h.c,
\end{equation}
where $g_{ik}$ is the strength of coupling between the $i$th level and $k$th mode. The coupling between the system and cavity mode is
\begin{equation}
\hat{V}_{sc}^{}=g(\hat{a}_l^\dag\hat{B}_{ba}+\hat{B}_{ba}^\dag\hat{a}_l).
\end{equation}
The strength of coupling between the system and the unimodal cavity
is denoted by $g$.
The system absorbs heat in the form of quanta from the hot bath and releases it in the cold bath after undergoing a radiative transition from $|a\rangle\rightarrow|b\rangle$.
Emission of coherent photons into the cavity as a result of this transition is the work done.

\subsection{Quantum Markovian Master Equation}
\label{s2a}
Within an adiabatic approximation, a weakly coupled quantum system's dynamics is governed by the quantum Liouville equation,
$
 |\dot\rho(t)\rangle=\breve{\mathcal L}(t)|\rho(t)\rangle,
$
where $|\dot\rho(t)\rangle$ is the time rate of change of the reduced density vector for the
system\cite{njp}.
$\breve{\mathcal L}(t)$ is the Liouvillian superoperator containing the time dependent driving and is responsible for system's evolution.
Using standard perturbation theory within the Born-Markov approximation,
the reduced density vector for the QHE is $|\rho\rangle=\{\rho_{11},\rho_{22},\rho_{aa},\rho_{bb},\Re (\rho_{12})\}$, where $\rho_{ii},i=1,2,a,b$ represent populations
 of the system's many body states and $\Re(\rho_{12})$ is the thermally (noise) induced coherence between the degenerate states $|1\rangle$ and $|2\rangle$.
   We assume that  the bath relaxation is much faster as compared
    to the driving time.
    In order to quantify the flux into the cavity mode, we
    focus on the statistics of the number of photons exchanged between the system and the  cavity.
    Within a measurement window, $t$,
the PDF corresponding to $q$ net photons in the cavity is $P(q,t)$. The statistics of $q$ is quantified by the
 moment generating function, defined as $G(\lambda,t)=\sum_qe^{\lambda q}P(q,t)$. From a general framework of FCS, it can be shown that
 the equation of motion for $G(\lambda,t)$ is
$
   \dot G(\lambda,t)=\langle \breve{{\boldsymbol 1}}| \breve{\mathcal L}(\lambda,t)|\rho(\lambda,t)\rangle
$\cite{dhar2,Levitov,uhrmp,uhFTfcs}.
 $\lambda$ is the auxiliary field that counts the number of photons exchanged between the system and the cavity.
  The transformed characteristic counting
  Liouvillian is given by (appendix), $\breve{\mathcal L}(\lambda,t)=$

  \begin{equation}
  \label{Liouvillian-lambda}
r\begin{pmatrix}
n(t)&0&\tilde{n}_h(t)&\tilde{n}_c(t)&y(t)\\
0&n(t)&\tilde{n}_h(t)&\tilde{n}_c(t)&y(t)\\
n_h(t)&n_h(t)&\frac{-g^2\tilde{n_l}-2r\tilde{n}_h^{}(t)}{r}&\frac{g^2n_l e^{-\lambda}}{r}&2p_hn_h(t)\\
n_c(t)&n_c(t)&\frac{g^2\tilde{n_l}e^{\lambda}}{r}&\frac{-g^2n_l-2r\tilde{n}_c(t)}{r}&2p_cn_c(t)\\
\frac{y(t)}{2}&\frac{y(t)}{2}&p_h\tilde{n}_h(t)&p_c\tilde{n}_c(t)&n(t)-\tau
\end{pmatrix}.\\[2mm]
\end{equation}
We have  denoted the time and counting-field dependent density vector as
$|\rho(\lambda,t)\rangle$ and $\langle\breve{{\boldsymbol 1}}|=\{1,1,1,1,0\}$. All
couplings between the QHE and thermal baths are
assumed to be equal and denoted by $r$. The adiabatic limit is valid for $rt'\gg 1$, where
$t'$ is the time scale of the external drivings.  Also, $\breve {\mathcal L}(\lambda=0,t)=\breve{\mathcal L}(t).$
Here, $n(t)=-n_c(t)-n_h(t)$, $y(t)=-n_c(t)p_c-n_h(t)p_h$, with
  $n_x(t),x\in h,c$, $\tilde{n}_x(t)=n_x(t)+1$. These are given by,
$
 n_c(t)=(\exp\{(E_b-E_1)/k_BT_c(t)\}-1)^{-1}, n_h(t)=(\exp\{(E_a-E_1)/k_BT_h(t)\}-1)^{-1}
$.
 The occupation of the cavity mode is
$
 n_l=1/(\exp[(E_a-E_b)/k_BT_l]-1), \tilde n_l=1+n_l
$, $T_l$ being a fictitious temperature of the cavity\cite{goswami2013thermodynamics}.
The dimensionless parameters, $p_h$ and $p_c$ represent quantum coherence control parameters associated with the hot and cold
baths respectively\cite{goswami2013thermodynamics,SaarQHEPhysRevA.86.043843,ScullyPNAS13092011}. These two parameters are a measure of interference
between the transitions involving the ground states and the two excited states.
$\tau$ is a dimensionless dephasing rate introduced phenomenologically \cite{ScullyPNAS13092011,SaarQHEPhysRevA.86.043843} so as to take care of environmental dephasing effects.

\subsection{Pancharatnam-Berry Curvature}
\label{s2b}
In the long time limit, PBp contributions can be realized in
the scaled cumulant generating function given by $S(\lambda)=\lim_{t\rightarrow\infty}(1/t)\ln G(\lambda,t)$.
 $S(\lambda,t)$ is
 additively separable into two parts\cite{sin} (appendix)
  $S(\lambda,t)=S_d(\lambda,t)+S_g(\lambda,t)$,
 \begin{eqnarray}
 \label{s-dyn}
S_d(\lambda)&=&\frac{1}{t_p}\displaystyle\int_0^{t_p}dt'\zeta_o(\lambda,t'),\\
S_g(\lambda)&=&-\frac{1}{t_p}\int_0^{t_p}\langle L_o({\lambda,t})
 |\dot R_o({\lambda,t})\rangle dt.
 \label{s-geo1}
 \end{eqnarray}
 Here, $S_{d(g)}(\lambda)$ denotes the dynamic(geometric) cumulant generating function. $|R_o(\lambda,t)\rangle $ and $\langle L_o(\lambda,t)|$ denote the instantaneous
 right and left eigenvectors of $\breve{\cal L}(\lambda,t)$ corresponding to the instantaneous long-time dominating
 eigenvalue, $\zeta_o(\lambda,t)$. Converting the line integral to a contour integral in the parameter space of $T_c(t)$ and $T_h(t)$ over a contour ${\cal C}$ and
 assuming
 the contour to be closed (a
 fixed time period) and piecewise smooth, we
can use Stokes' theorem and
rewrite the contour integral as a surface integral
 enclosed by the closed loop $\mathcal S$,
 \begin{align}
\label{s-geo}
 S_g(\lambda)&=-\frac{1}{t_p}\oiint_{\mathcal S}\nabla\times\langle L_o({\lambda, \bf T})
 |\partial_{{\bf T}}|R_o({\lambda,
 \bf T})\rangle. d{\mathcal S}.\nonumber\\
 \end{align}
 Here, vector ${\bf T}$ contains system parameters modulated by the external driving in the contour.
Eq. (\ref{s-geo})  has a
geometric interpretation since it is quantified by the surface's area and is both re-parametrization as well as gauge invariant\cite{sin, ren2010berry,hpg4}.
The geometric interpretation of Eq. (\ref{s-geo}) is analogous to the original GP interpretation in
 isolated quantum dynamics \cite{berry1984quantal,mukunda1993quantum,aharonov1987phase}, albeit not being a phase.
 We simply refer to it as the geometric phase-like or Pancharatnam-Berry phase-like  contribution.
The measurement time,
 $t=n t_p$, where $n$ is the number of cycles and
  $t_p$ is the time-period of the driving such that $r t_p\gg 1$.
$\langle L_o({\lambda, \bf T})
 |\partial_{{\bf T}}|R_o({\lambda,
 \bf T})\rangle$ is equivalent to
 the geometric potential in the parametric  space of $T_c(t)$ and $T_h(t)$ whose curl gives the curvature. It vanishes when there is no phase-difference between
 the driving protocols (also for $\phi$ as integer multiples
  of $\pi$) since there is no area traced
  in the parameter space.

\section{Thermodynamics}
\label{s3}

      \begin{figure}
\centering
\includegraphics[width=8.5cm]{fig_1.pdf}
\caption{  a) Dynamic flux contour. Dashed line indicates that $p_{h,d}^*$ is constant along $\phi$.
b) The geometric flux contour.
The optimal values occur at zero coherences. c) $j/j_o$ contour (dimensionless). Due to higher geometric contribution, the optimal flux values
occur at zero coherence value. $E_1=E_2=0.1,E_b=0.4,E_a=1.5,r=0.7,g=10,\tau=0.5,p_c=0.4,\omega=1.2,T_h=4,T_c=2,T_l=1,A=0.2.$
d) Linear (no) dependence of geometric (dynamic) flux on driving frequency.
 Energies ($E_1,E_2,E_b,E_a$), temperatures ($T_h,T_c,T_l$) are in atomic units ($\hbar=1,k_B=1$), rates $g^2,r$ have the
 units of time inverse and the driving frequency, $\omega$ is in $kHz$ throughout all figures.
}
\label{jdyn}
\end{figure}
We now focus on the flux, work and efficiency of the QHE.
The flux is defined as the rate of change of the number of photons exchanged between the system and the unimodal cavity. It can be obtained
 using, $j=\partial_\lambda S(\lambda)|_{\lambda=0}$\cite{uhrmp}. Since $S(\lambda)$ is separable into a geometric and dynamic part,
  the flux is a sum of   dynamic ($j_d$) and geometric ($j_g$) contributions.
  For sinusoidal drivings $T_h(t,\phi)=T_h(1-A\sin(\omega t+\phi))$ and $T_c(t)=T_c(1-A\sin(\omega t))$, with the amplitude $A<1$, $\omega$ as the driving frequency and
   $\phi$ as the phase difference between the two driving protocols, we
   numerically  evaluated the flux which is shown in Fig.(\ref{jdyn}).
 Similar to a non-driven engine\cite{goswami2013thermodynamics}, the dynamic
 flux can be maximized with respect to the coherence parameter $p_h$ (dashed line in Fig.(\ref{jdyn}a) for a fixed $p_c=0.4$).
   Note that, this coherence-maximized value of the dynamic flux (at constant $p_c$) as a function of phase difference is always at a constant $p_h$ value, i.e at $p_h=p_{h,d}^*$.
  In Fig. (\ref{jdyn}a), at $p_h=p_{h,d}^*(=0.35)$, the coherence-maximized value of the flux decreases and reaches a minimum value at $\phi=\pi$
  (where the two drivings are completely out of phase and the thermodynamic force is minimum). Infact for a fixed value of $p_h$, the flux always decreases as a function of $\phi$,
  being minimum at $\phi=\pi$ before increasing again. This trend is robust for all values of $p_c$ and $p_h$.
  Also,
    $p_{h,d}^*$ for a fixed set of other parameters
     linearly depends on $p_c$ and this linearity is also independent of the phase difference as shown
    in Fig.(\ref{ph-linear}a).  This linear dependence of $p_{h,d}^*$ on $p_c$ was also observed in the nondriven heat engine \cite{goswami2013thermodynamics}.

    \begin{figure}
\centering
\includegraphics[width=8.5cm]{fig_2.pdf}
\caption{a) Plot showing the linear (no) dependence of $p_{h,d}^* (p_{h,g}^*)$ on $p_c$ and both
 are independent of the phase difference ($p_{h,g}^*$ is the dashed magenta line).
b) Near equilibrium EMP as a function of Carnot efficiency. Only for $\phi=0,\pi$, the slope is 1/2.
For other phase difference, linear slopes greater ($\phi=\pi/3$, slope=0.65) and lower ($\phi=4\pi/3$, slope=0.42)
than 1/2 are observed. Dashed black (cyan) line is the lower (upper) bounds on EMP, $T_c=0.9$, $T_l=2$, $p_c=0.1$, $p_h=1$, $g=50$, $\omega=3,$ $A=0.02, r=0.7,\tau=0.5$.
}
\label{ph-linear}
\end{figure}

  The interesting quantity, however, is the geometric flux.
  It doesn't typically follow the same dependence on the phase difference as the dynamic flux does. As can be seen from the contour Fig. (\ref{jdyn}b),
   the geometric flux, oscillates as a function of the phase difference. For any arbitrary $p_h$, we see that it increases from zero ($\phi=0$
   ) and reaches a maxima at $\phi=\pi/2$ before decreasing to become zero again at $\phi=\pi$. After that it is negative, reaches a minima at $\phi=3\pi/2$ and
   then increases to zero at $\phi=2\pi$.
   The negative component of the geometric flux has been utilized in creating pumping in quantum transport \cite{ren2010berry,hayakawa}.
   In our case, we further observe that the value of $p_h$ that optimizes  the geometric flux, $p_{h,g}^*$,
   is always at $p_{h,g}^*=0~\forall ~p_c,\phi$ as
    shown in Fig. (\ref{ph-linear}a, dashed magenta line).
      This  indicates
     that coherences are exclusive to optimizing the  dynamic flux alone.
  Since the geometric flux increases linearly with $\omega$, Fig.(\ref{jdyn}d) and the dynamic flux
      doesn't depend on $\omega$, the contribution to total flux is controllable through $\phi$ and $\omega$.
           In the limit of high $\omega$ and nonzero $\phi (\phi\ne n\pi,n\in$ integers), $j_g$  contributes dominantly to the total flux and
     the global maxima (minima) of the total flux occurs when there is no coherence. This is shown in Fig. (\ref{jdyn}c), where the global maxima (minima)
     of $j/j_o$ ($j_o$ is the flux  when $\phi=0, p_c=p_h=0$ or $ p_c=p_h$)
      is at $p_h=0,\phi=\pi/2(3\pi/2)$. In this limit, the optimum value of the total flux occurs in absence
       of quantum coherences.
     Quantum coherences are known to optimize the flux and power of QHEs beyond classical values\cite{ScullyPNAS13092011} in quantum heat engines.
      However, by increasing the geometric flux, one can nullify the effect of quantum coherences in optimizing the total flux by increasing the frequency of driving.
      Note that the upper limit of $\omega$ is in the $THz$ regime beyond which the adiabatic approximation will fail\cite{ren2010berry}.
 Further, changes in the decoherence parameter, $\tau$, only cause slight magnitude shifts in the flux and doesn't change the inference of the observations.

 In the cavity, coherent photons of energy $E_{ab}=E_a-E_b$  are generated each time the system relaxes from state $|a\rangle$ to $|b\rangle$.
 However,
there is dissipation in the cavity mode due to stimulated emission. Since the cavity occupation, $n_l$, has to be kept constant,
this dissipation is proportional  to $\ln(\tilde n_l/n_l)$
\cite{goswami2013thermodynamics}. By taking care of this dissipation, the actual work done by the QHE can be written as,
\begin{eqnarray}
 W= E_{ab}-\alpha\ln\frac{\tilde{n}_l}{n_l},
\end{eqnarray}
with $\alpha$ as the proportionality factor. For a QHE with no driving $\alpha=k_BT_c$ \cite{goswami2013thermodynamics}.
 In this case, since there is a cyclic driving, the proportionality factor will depend on the driving protocol. However, per cycle,
  the proportionality factor will be $\alpha=k_B/t_p\int_0^{t_p}T_c(t)dt$.
  This is because, the cavity dynamics is equivalent for both the driven and non-driven cases with the only difference $T_c\rightarrow T_c(t)$. Using this definition,
 we calculate the efficiency at maximum power (EMP) and focus on the expansion coefficients near equilibrium.
 In order to calculate the EMP we maximize the power, $P=jW$ with respect to the energy $E_b$.
 Near equilibrium, we observe that the EMP goes linearly, $\eta_*=m\eta_c$.
 It has been shown that $m=1/2$  and this value is universal in the close to equilibrium in the endoreversible regime\cite{PhysRevLett.102.130602,BroeckEMP}.
 However, we see that $m=1/2$ only for $\phi$ being integer multiples of $\pi$ (i.e no PBp effects).
 For other phase differences $\phi>(<)\pi$, $m<(>)1/2$  as shown in Fig. (\ref{ph-linear}b). We further observe that the EMP is not restrictive
  in the bounds $\eta_c/2\le \eta_*\le\eta_c/(2-\eta_c)$\cite{izumida2012efficiency,PhysRevE.87.012133} as seen from Fig. (\ref{ph-linear}b). When $\phi< (>)\pi$, the EMP
  crosses the upper (lower)  bound on $\eta_*$.
  This can be explained from the oscillatory behavior of the geometric flux. Since $j_g>(<)0 \forall \phi<(>)\pi$, it adds (subtracts) to (from) the total flux, thereby
   changing the position of $E_b$ where $P$ is maximum,  increasing (decreasing) the slope of $\eta_*$.
 Deviations from universality in near equilibrium expression of $\eta_*$ has also been recently discussed in a classical heat engine with heat exchangers in the linear response regime\cite{PhysRevE.96.022119}.

\section{Steady State Fluctuation Theorem}
\label{s4}
In a general framework of open quantum systems,
the FTs \cite{dhar2,uhrmp,campisi2011colloquium,utsumi2010bidirectional} have its manifestation in the generating function, $G(\lambda,t)$.
   A symmetry of the type $G(\lambda,t)=G(-\lambda- F,t)$, popularly called the Gallavoti-Cohen (GC) symmetry \cite{Galla-Cohen,dhar2,uhrmp,max-uh} exists
   and is analogous to the FT. Any open quantum system with the GC symmetry satisfies this standard form of the SSFT \cite{uhrmp,max-uh}.
In open quantum systems, SSFT in the  long observation time reads,
$
 \ln\lim_{t\rightarrow\infty}[P(q,t)/P(-q,t)]=qF,
$
where $F$ is the thermodynamic
 force that drives the system out of equilibrium and has recently been experimentally verified in quantum dots \cite{utsumi2010bidirectional}.
 The FT guarantees positive entropy production since the probability of observing negative flux is exponentially suppressed.
  The PDF is connected to the long time scaled generating function by the Gartner-Ellis-Varadhan theorem
 \cite{touchette2009large,varadhan1984large},
$
  P(q,t)\asymp \exp\{-t L(y)\}$ with $
  L(y)=\sup_\lambda^{}(y\lambda-S(\lambda)),$
 $L(y)$ being the large-deviation function (LDF), a Legendre transform of the scaled cumulant generating function with  $y=q/t$. The LDF obeys the symmetry,
$
 L(y)-L(-y)=-yF.
$ The symmetry in the LDF and the GC symmetry are implications of the SSFT\cite{uhrmp,dhar2}.

    \begin{figure}
\centering
\includegraphics[width=8.5cm]{fig_3.pdf}
\caption{a) GC symmetry  $S_d(\lambda) $(red)$=S_d(-\lambda-F)$ (dotted), $F=0.235$. b) Symmetry in the dynamic LDF with a slope of $-F$. c) Asymmetric behavior
 of $S_g(\lambda)$ which  breaks SSFT. d) Non-convex  nature of the linearly shifted $S_g(\lambda)$.
 Parameters: $T_c=1,T_h=2.5,T_l=2,A=0.15,\omega=1.2, p_c=p_h=0,\phi=\pi/3,r=0.7,g=50,\tau=0.5$.}
\label{GCfig}
\end{figure}

In the QHE, we evaluate the GC symmetry by numerically calculating the
 eigenvectors of the counting Liouvillian, Eq. (\ref{Liouvillian-lambda}). We find that the GC symmetry
 holds only for the dynamic part, $S_d(\lambda)=S_d(-\lambda-F)$, Fig. (\ref{GCfig}a) as seen in some previous works
  with both heat and electron transport
 \cite{ren2010berry,hpg4} for any phase difference. Likewise, the symmetry in the LDF holds only for the dynamic part $L_d(y)-L_d(-y)=-yF$, Fig. (\ref{GCfig}b) .
 The thermodynamic force in the QHE is (appendix),
 \begin{eqnarray}
 F=\displaystyle\ln \frac{\tilde n_l
 \int_0^{t_p}\tilde n_c n_hdt}{n_l\int_0^{t_p}n_c\tilde n_hdt}.
 \end{eqnarray}
 However, for a nonzero phase difference ($\phi\ne \pi,2\pi\ldots$),  $S_g(\lambda)$ also contributes to the PDF.
Unlike the dynamic part, the GC symmetry doesn't hold for the geometric part, $S_g(\lambda)\ne S_g(-\lambda-R)$ (Fig. \ref{GCfig}c). Infact, there is no symmetry in $S_g(\lambda).$
This causes  the form of the SSFT to be violated, a direct consequence of finite PBp effects.
Further,
  $S_g(\lambda)$ as well as the linearly shifted function, $y\lambda - S_g(\lambda)$, (Fig. \ref{GCfig}d), are not  strict monotonously increasing functions.
Hence, the standard Gartner-Ellis-Varadhan theorem cannot be used to evaluate the PDF anymore since a standard Legendre transformation cannot be performed. This inapplicability is due to the fact that the
LDF is based on a saddle-point approximation valid for convex upwards or concave downward functions \cite{touchette2009large, varadhan1984large},
 which is  not the case here.
Note that, the PDF  can be obtained by doing a standard inverse transform on the GF.
We conclude that the standard form of the SSFT and the applicability of the large deviation theory is restricted to the case when
PBp contributions are absent, i.e the curvature in Eq. (\ref{s-geo}) is zero.

\section{Conclusions}
\label{s5}
We investigated a periodically driven quantum heat engine by parametrically modulating the temperatures of the reservoir
within an adiabatic quantum Markovian master equation approach.
From a general framework of full counting statistics, we showed that the scaled cumulant generating function corresponding to the number of photons emitted in the cavity
is separable into  dynamic and Pancharatnam-Berry phase-like
contributions.
The Pancharatnam-Berry phase-like contributions were shown to violate the steady state fluctuation theorem as well as the expansion coefficient and the  bounds on the universal expression
 for the efficiency at maximum power. We also showed that one cannot use
  a large deviation technique to evaluate the probability distribution function.
  These violations can however be recovered for a vanishing geometric curvature. Further we observed that,
  the effect of quantum coherences in optimizing the total flux can be nullified by having a dominant geometric contribution to the total flux.

\begin{acknowledgments}
   We acknowledge the support from the Max-Planck-Institute for the Physics of Complex Systems, Dresden, Germany and thank
   Izaak Neri for going through the manuscript.
   HPG would like to thank Prof. Udo Seifert for a helpful discussion.
 \end{acknowledgments}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{sec:introduction}
In the framework of Color Glass Condensate (CGC) effective field theory (EFT)~\cite{Weigert2005,Jalilian-Marian2005,Iancu2003,Gelis2010}, both projectile and target in high energy collisions are regarded as high parton density sources. The effective degrees of freedom are color sources $\rho$ at large $x$ and classical gauge fields $\mathcal{A}_\mu$ at small $x$. The classical gauge field $\mathcal{A}_\mu$ is the solution of classical Yang-Mills equations with a fixed configuration of color sources at an initial scale $x_0$. Its dynamical evolution is captured by the Balitsky-Kovchegov (rcBK) equation with running coupling kernel~\cite{Balitsky1996,Balitsky1999,Kovchegov1999}. Gluon distribution in this regime is saturated with typical transverse momentum $\vert\bo{k}_\perp\vert\sim Q_s$, and localized in typical size $\sim 1/Q_s$. $Q_s$ is the saturation scale.

A strongly interacting matter (Glasma) is formed after a short time of collisions by the two incident particles with CGC dynamics.  Multi-particle production is a consequence of an approximate boost invariant radiation from Glasma flux tubes. Using the CGC EFT and combining the fragmentation models of hadronization, the multiplicity distributions in $pp$ and $pA$ collisions at Relativistic Heavy Ion Collider (RHIC) through Large Hadron Collider (LHC) energies are successfully described~\cite{Dumitru2006a,Dumitru2006b,Goncalves2006,Tuchin2008,Albacete2010,Rezaeian2010,
Fujii2011,Albacete2012,Chirilli2012,Chirilli2012a}. The transverse momentum distributions of charged particles and broadening predicted by the CGC as a function of multiplicity is also clearly seen in the data~\cite{Albacete2012,McLerran2010}.

In particular, the data of dihadron correlations at RHIC and LHC~\cite{Adams2005,PHENIXCollaboration2008,Alver2009,CMSCollaboration2010,CMSCollaboration2012a,
ALICECollaboration2012,CMSCollaboration2013} are qualitatively understood by some of the phenomenological models with the early-time dynamics determined by CGC EFT\cite{Dumitru2008c,Dusling2009,Dumitru2011d,Gelis2008,Dumitru2010a,
Dusling2012b,Dusling2013d,Dusling2013e,Dusling2013b}. In our previous paper~\cite{Zhao2015}, we calculate the two-gluon rapidity correlations of $pp, pA$ and $AA$ collisions in the framework of CGC EFT to leading order of $p_\perp/Q_s$ under fixed coupling constant $\alpha_s$. The results are qualitatively consistent with current data for $pPb$ collisions presented by the ALICE Collaboration~\cite{Alice2016}. It is firstly noted that the ridge-like long-range rapidity correlations are caused by the stronger correlation with the gluon of colour source. The ridge is more likely observed at higher incident energy and lower transverse momentum of trigger gluon.

To see quantitatively the influences of the sub-leading order and running couple to two-gluon rapidity correlations, in this paper, we will calculate the contributions of the leading and the sub-leading orders of $p_\perp/Q_s$ under running coupling constant $\alpha_s$. It will present a more complete and accurate description of CGC EFT on two-gluon rapidity correlation pattern.

The paper is organized as follows. The formulas of single and double-gluon inclusive production amplitudes in the leading and sub-leading orders are firstly clarified in section~\ref{sec:review}.  Then numerical results of two-gluon rapidity correlations in the leading and the sub-leading orders are presented and discussed respectively in section~\ref{sec:results}. Finally, a brief summary and conclusions are given in section~\ref{sec:summary}.

\section{Two-gluon correlator in the framework of CGC}
\label{sec:review}
Supposing two gluon produced with transverse momentum
$\bo{p}_\perp$ and $\bo{q}_\perp$, and the longitudinal rapidity
$y_p$ and $y_q$. Their correlator is usually defined as,
\begin{equation}\label{eq:24}
C(\bo{p}_\perp,y_p;\bo{q}_\perp, y_q)= \frac{
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp
dy_q}\right\rangle }{ \left\langle\frac{dN_1}{d^2\bo{p}_\perp
dy_p}\right\rangle \left\langle\frac{dN_1}{d^2\bo{q}_\perp
dy_q}\right\rangle }-1=\frac{
\left\langle\frac{dN^{\mathrm{corr.}}_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle}{\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle
\left\langle\frac{dN_1}{d^2\bo{q}_\perp dy_q}\right\rangle},
\end{equation}
where $\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp
dy_q}\right\rangle$ and $\left\langle\frac{dN_1}{d^2\bo{p}_\perp
dy_p}\right\rangle$ are the double and single gluon inclusive
productions.
$\left\langle\frac{dN^{\mathrm{corr.}}_2}{d^2\bo{p}_\perp dy_p
d^2\bo{q}_\perp dy_q}\right\rangle$ is the two correlated gluon production, where the uncorrelated part is subtracted.

In the framework of CGC EFT~\cite{Gelis2008}, for a given collision, the observable under the leading log approximation is factorized as~\cite{Gelis2008},
\begin{equation}\label{eq:llog}
\langle\mathcal{O}\rangle_{\mathrm{LLog}}=\int\left[D\rho_1\right]\left[D\rho_2\right]W[\rho_1]W[\rho_2]\mathcal{O}[\rho_1,\rho_2]_{\mathrm{LO}}.
\end{equation}
Here $\mathcal{O}[\rho_1,\rho_2]_{\mathrm{LO}}$ is the leading order single or double gluon inclusive distribution for a fixed distribution of color sources, and the integration denotes an average over different distribution of the color sources with the weight functional $W[\rho_{1,2}]$. In general, $W[\rho_{1,2}]$ encodes all possible color charge configurations of projectile and target, and obeys Jalilian-Marian-Iancu-McLerran-Weigert-Kovner (JIMWLK) renormalization group equations~\cite{Balitsky1996,Jalilian-Marian1997,Jalilian-Marian1998,Kovner2000,Iancu2001,Iancu2001a,Ferreiro2002}.
 Where all quantum information of projectile/target is absorbed into the distribution $W[\rho_{1,2}]$.

Let's firstly take the simplest case, the single gluon inclusive production, to show the contribution of leading order~\cite{Blaizot2004b},
\begin{equation}
\left.\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\vert_{\mathrm{LO}}=\frac{1}{2(2\pi)^3}\sum_{a,\lambda}\lvert\mathcal{M}^a_\lambda(\bo{p})\rvert^2,
\end{equation}
where the amplitude can be written as
\begin{equation}\label{eq:single-amplitude}
\mathcal{M}^a_\lambda (\bo{p})= p^2A^{a,\mu}(\bo{p})\epsilon^{(\lambda)}_\mu (\bo{p}),
\end{equation}
here, $p$ and $\epsilon^{(\lambda)}_{\mu}(\bo{p})$ are the 4-momenta and polarization vector of produced gluon, respectively. $A^{a,\mu}$ with color index $a$ is the background gauge field, which is the solution of the classical Yang-Mills equation with a fixed configuration of color charge sources for projectile/target. The amplitude for this process is shown in~\fig{fig:single-gluon}.
\begin{figure}[hbt]
\begin{center}
\includegraphics[scale=0.8]{figs/single.eps}
\end{center}
\caption{The diagram of the single-gluon inclusive production in a collision. The vertical dashed line denotes final state.}
\label{fig:single-gluon}
\end{figure}

Performing an average over all possible configurations of color charges sources of projectile/target in \eq{eq:llog}, we get,
\begin{equation}\label{eq:single-llog}
\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle_{\mathrm{LLog}}=\frac{1}{2(2\pi)^3}\sum_{a,\lambda}\left\langle\lvert\mathcal{M}^a_\lambda(\bo{p})\rvert^2\right\rangle_{A_{1,2}}.
\end{equation}
At large transverse momentum $k_\perp\gg Q_s$ ($Q_s$ is the saturated scale of incident nucleus), the background gauge field can be expanded in power of $\tilde{\rho}/\bo{k}^2_\perp$\cite{Blaizot2004b,Dusling2009,Kovner1995a,Kovchegov1997},
\begin{equation}\label{eq:expand}
p^2 A^{a,\mu}(\bo{p})= -if^{abc}g^3\int\frac{d^2\bo{k}_\perp}{(2\pi)^2}L^\mu(\bo{p},\bo{k}_\perp)\frac{\tilde{\rho}^b_1(\bo{k}_\perp)\tilde{\rho}^c_2(\bo{p}_\perp-k_\perp)}{\bo{k}^2_\perp(\bo{p}_\perp-\bo{k}_\perp)^2},
\end{equation}
Where $k_\perp=\vert\bo{k}_\perp\vert$, $f^{abc}$ is the structure constant of $SU(3)$, $\tilde{\rho}_{1,2}$ are the Fourier modes of color charge density of projectile and target, respectively. $L^\mu(\bo{p},\bo{k}_\perp)$ is the Lipatov vertex\cite{Kuraev1977,Balitsky1978}. In light-cone coordinate, it is expressed explicitly as
\begin{eqnarray}
\begin{split}
L^+(\bo{p},\bo{k}_\perp) &= -\bo{k}^2_\perp /p^-,\\
L^-(\bo{p},\bo{k}_\perp) &=
\left[(\bo{p}_\perp-\bo{k}_\perp)^2-\bo{p}^2_\perp \right]/p^+,\\
L^i(\bo{p},\bo{k}_\perp)&=-2\bo{k}^i_\perp.
\end{split}
\end{eqnarray}
Substituting \eq{eq:single-amplitude} and \eq{eq:expand} into \eq{eq:single-llog}, and summing over the polarizations and color indices, the single-gluon inclusive production in LLog approximation becomes,
\begin{equation}\label{eq:7}
\begin{split}
\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle_{\mathrm{LLog}}=&-\frac{1}{2(2\pi)^3}f^{abc}f^{ade}g^6\int\frac{d^2\bo{k}_\perp}{(2\pi)^2}\frac{d^2\bo{k}'_\perp}{(2\pi)^2}L^\mu (\bo{p},\bo{k}_\perp)L_\mu (\bo{p},\bo{k}'_\perp)\\
&\times
\frac{\left\langle
\tilde{\rho}^b_1(y_p,\bo{k}_\perp)
\tilde{\rho}^c_2(y_p,\bo{p}_\perp-\bo{k}_\perp)
\tilde{\rho}^{\ast d}_1(y_p,\bo{k}'_\perp)
\tilde{\rho}^{\ast e}_2(y_p,\bo{p}_\perp-\bo{k}'_\perp)\right\rangle_{A_{1,2}}
}{\bo{k}^2_\perp (\bo{p}_\perp-\bo{k}_\perp)^2
\bo{k}'^{2}_\perp (\bo{p}_\perp-\bo{k}'_\perp)^2}.
\end{split}
\end{equation}
The average $\langle\cdots\rangle$ in \eq{eq:7} is over all possible color charge configurations of projectile and target. In the McLerran-Venugopalan (MV) model\cite{McLerran1994,McLerran1994a,McLerran1994b}, the color charge distribution in each nucleus has a local Gaussian form, the average of multi color charge densities in \eq{eq:7} can be expressed in terms of the two-point correlator,
\begin{equation}\label{eq:8}
\langle\tilde{\rho}^a(\bo{k}_\perp)\tilde{\rho}^{\ast b}(\bo{k}'_\perp)\rangle_{A_{1,2}}=(2\pi)^2\delta^{ab}\delta^2(\bo{k}_\perp -\bo{k}'_\perp)\mu^2_{A_{1,2}}(y_p,\bo{k}_\perp),
\end{equation}
here $\mu^2(y_p,\bo{k}_\perp)$ is the color charge squared per unit transverse area in momentum space, $y_p$ is the longitudinal rapidity of produced gluon. The \eq{eq:7} becomes,
\begin{equation}\label{eq:9}
\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle_{\mathrm{LLog}}=\frac{-g^6S_\perp N_c(N^2_c-1)}{2(2\pi)^3}\int\frac{d^2\bo{k}_\perp}{(2\pi)^2}
\frac{L^\mu(\bo{p},\bo{k}_\perp)L_\mu(\bo{p},\bo{k}_\perp)}{\bo{k}^4_\perp (\bo{p}_\perp-\bo{k}_\perp)^4}
\mu^2_{A_1}(y_p,\bo{k}_\perp)\mu^2_{A_2}(y_p,\bo{p}_\perp - \bo{k}_\perp).
\end{equation}
Where $f^{abc}f^{ade}=N_c(N^2_c-1)$. The un-integrated gluon distribution (UGD)~\cite{Dumitru2008c,Dusling2009} is given by,
\begin{equation}\label{eq:10}
\Phi_{A_{1,2}}(y_p,\bo{p}_\perp)=g^2\pi(N^2_c-1)\frac{\mu^2_{A_{1,2}}(y_p,\bo{p}_\perp)}{\bo{p}^2_\perp}.
\end{equation}
Substituting \eq{eq:10} into \eq{eq:9} and following some algebras, the \eq{eq:9} becomes,
\begin{equation}\label{eq:11}
\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle_{\mathrm{LLog}}=\frac{\alpha_s N_c S_\perp}{\pi^4(N^2_c-1)}\frac{1}{\bo{p}^2_\perp}\int\frac{d^2\bo{k}_\perp}{(2\pi)^2}\Phi_{A_1}(y_p,\bo{k}_\perp)\Phi_{A_2}(y_p,\bo{p}_\perp - \bo{k}_\perp).
\end{equation}
\eq{eq:11} is the final result of the single-gluon inclusive distribution applied LLog factorization approach, which is also well-known $k_\perp-$factorization formalism~\cite{Kovchegov1998} and widely used in the computation of single inclusive particle production at high energy collisions. $S_\perp$ is the transverse overlap area of two UGDs, i.e., the integral of two $\delta$-functions. And the running coupling constant $\alpha_s$ in \eq{eq:10} and \eq{eq:11} changes with $p_\perp$.  It is fixed in our previous paper~\cite{Zhao2015}.

Next, multi-gluon inclusive distribution is straightforward to derive in this way. For specific, double-gluon inclusive production, if we label the two produced gluon as $p$ and $q$, transverse momentum $\bo{p}_\perp$ and $\bo{q}_\perp$, and the longitudinal rapidity $y_p$ and $y_q$, the double-gluon distribution in the LLog approximation can be directly written as,
\begin{equation}\label{eq:12}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=\frac{1}{2^2(2\pi)^6}
\sum_{a,a';\lambda,\lambda'}
\left\langle\lvert\mathcal{M}^{aa'}_{\lambda\lambda'}(\bo{p},\bo{q})\rvert^2\right\rangle_{A_{1,2}},
\end{equation}
where the amplitude is,
\begin{equation}\label{eq:13}
\mathcal{M}^{aa'}_{\lambda\lambda'}(\bo{p},\bo{q})
=p^2q^2A^{a,\mu}(\bo{p})A^{a',\sigma}(\bo{q})\epsilon^{(\lambda)}_{\mu}(\bo{p})\epsilon^{(\lambda')}_{\sigma}(\bo{q}).
\end{equation}
$\epsilon^{(\lambda)}_{\mu}(\bo{p})$ and $\epsilon^{(\lambda')}_{\sigma}$ are the polarization vectors of two produced gluon. Substituting \eq{eq:expand} and \eq{eq:13} into \eq{eq:12}, one arrives at,
\begin{equation}\label{eq:14}
\begin{split}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=&
\frac{1}{4(2\pi)^6}
f^{abc}f^{a'de}f^{a\bar{b}\bar{c}}f^{a'\bar{d}\bar{e}}
g^{12}\int\frac{d^2\bo{k}_{1\perp}}{(2\pi)^2}
\frac{d^2\bo{k}_{2\perp}}{(2\pi)^2}
\frac{d^2\bo{k}_{3\perp}}{(2\pi)^2}
\frac{d^2\bo{k}_{4\perp}}{(2\pi)^2}\\
&\times L^\mu(\bo{p},\bo{k}_{1\perp})
L^\sigma(\bo{q},\bo{k}_{3\perp})
L_\mu(\bo{p},\bo{k}_{2\perp})
L_\sigma(\bo{q},\bo{k}_{4\perp})\\
&\times
\frac{\left\langle
\tilde{\rho}^b_1(\bo{k}_{1\perp})
\tilde{\rho}^c_2(\bo{p}_\perp-\bo{k}_{1\perp})
\tilde{\rho}^{\ast \bar{b}}_1(\bo{k}_{2\perp})
\tilde{\rho}^{\ast \bar{c}}_2(\bo{p}_\perp-\bo{k}_{2\perp})
\tilde{\rho}^d_1(\bo{k}_{3\perp})
\tilde{\rho}^e_2(\bo{q}_\perp-\bo{k}_{3\perp})
\tilde{\rho}^{\ast\bar{d}}_1(\bo{k}_{4\perp})
\tilde{\rho}^{\ast\bar{e}}_2(\bo{q}_\perp-\bo{k}_{4\perp})
\right\rangle_{A_{1,2}}
}{
\bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2
\bo{k}^2_{2\perp}(\bo{p}_\perp-\bo{k}_{2\perp})^2
\bo{k}^2_{3\perp}(\bo{q}_\perp-\bo{k}_{3\perp})^2
\bo{k}^2_{4\perp}(\bo{q}_\perp-\bo{k}_{4\perp})^2
}.
\end{split}
\end{equation}
In \eq{eq:14}, one can evolute the quantities of the average over color sources by applying local Gaussian correlators, \textit{i.e.}, \eq{eq:8} for equal rapidity correlators. As for non-equal rapidity correlators, \textit{i.e.}, $y_q>y_p$
\begin{equation}\label{eq:15}
\begin{split}
\langle\tilde{\rho}^a_1(y_p,\bo{k}_\perp)\tilde{\rho}^b_1(y_q,\bo{k}'_\perp)\rangle
&=\langle\tilde{\rho}^a_1(y_p,\bo{k}_\perp)\left[\tilde{\rho}^b_1(y_p,\bo{k}'_\perp)+\Delta\tilde{\rho}^b_1(y_q-y_p,\bo{k}'_\perp)\right]\rangle\\
&=(2\pi)^2\delta^{ab}\delta^2(\bo{k}_\perp - \bo{k}'_\perp)
\mu^2_1(y_p,\bo{k}_\perp).
\end{split}
\end{equation}
Here we assume $y_q>y_p$. $\Delta\tilde{\rho}$ is the difference of color charge densities of two-gluon is ignored as at high energies, the color charge density inside each nucleus is saturated and the relative occupation number of gluons $\Delta y\alpha_s\ll 1$.

Now, $\langle\cdots\rangle_{A_{1,2}}$ can be evaluated by applying \eq{eq:8} and \eq{eq:15}. It contains 9 distinct diagrams. One is disconnected, and eight of them are connected. Let's demonstrate them one by one, respectively.
\begin{figure}[hbt]
\begin{center}
\includegraphics[scale=0.8]{figs/double-1-1.eps}
\end{center}
\caption{The disconnected diagram of the double-gluon inclusive production. The double solid lines corresponding to the transverse positions of scattered quarks in one of the incident nuclei.}
\label{fig:1-1}
\end{figure}
In \fig{fig:1-1}, we show the disconnected diagram for the double-gluon inclusive production. There is a spatial gap between the scattered quarks as denoted by double solid lines in \fig{fig:1-1}. Scattered quarks in each nuclei are localized in a transverse area of size $1/Q_s$, which is corresponding to that of \eq{eq:8}. Substituting \eq{eq:8} into \eq{eq:14}, one obtains,
\begin{equation}\label{eq:16}
\begin{split}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=&\frac{1}{4(2\pi)^6} f^{abc}f^{a'de}f^{a\bar{b}\bar{c}}f^{a'\bar{d}\bar{e}}g^{12}
\int\frac{d^2\bo{k}_{1\perp}}{(2\pi)^2} \frac{d^2\bo{k}_{2\perp}}{(2\pi)^2}
\frac{d^2\bo{k}_{3\perp}}{(2\pi)^2} \frac{d^2\bo{k}_{4\perp}}{(2\pi)^2}\\
&\times L^\mu(\bo{p},\bo{k}_{1\perp})L_\mu(\bo{p},\bo{k}_{1\perp})
L^\sigma(\bo{q},\bo{k}_{3\perp})L_\sigma(\bo{q},\bo{k}_{3\perp})\\
&\times\frac{(2\pi)^8\delta^{b\bar{b}}\delta^{d\bar{d}}\delta^{ce}\delta^{\bar{c}\bar{e}}
\mu^2_1(y_p,\bo{k}_{1\perp})\mu^2_1(y_q,\bo{k}_{3\perp})\mu^2_2(y_p,\bo{p}_\perp-\bo{k}_{1\perp})\mu^2_2(y_q,\bo{q}_\perp-\bo{k}_{3\perp})}
{\bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2 \bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2
\bo{k}^2_{3\perp}(\bo{q}_\perp-\bo{k}_{3\perp})^2
\bo{k}^2_{3\perp}(\bo{q}_\perp-\bo{k}_{3\perp})^2}.
\end{split}
\end{equation}
Following some algebras, \eq{eq:16} casts into
\begin{equation}\label{eq:17}
\begin{split}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=&\frac{\alpha^2_\mathrm{s} N^2_\mathrm{c}S_\perp}{\pi^8(N^2_\mathrm{c}-1)^2}\frac{1}{\bo{p}^2_\perp \bo{q}^2_\perp}\\
&\quad\times\int\frac{d^2\bo{k}_{1\perp}}{(2\pi)^2} \frac{d^2\bo{k}_{3\perp}}{(2\pi)^2}
\Phi_1(y_p,\bo{k}_{1\perp})\Phi_2(y_p,\bo{p}_\perp-\bo{k}_{1\perp})
\Phi_1(y_q,\bo{k}_{3\perp})\Phi_2(y_q,\bo{q}_\perp-\bo{k}_{3\perp})\\
=&\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle_\mathrm{LLog}
\left\langle\frac{dN_1}{d^2\bo{q}_\perp dy_q}\right\rangle_\mathrm{LLog}.
\end{split}
\end{equation}
It is the final result of the disconnected diagram as shown in \fig{fig:1-1}, which can be seen as the convolution of two single-gluon diagrams as shown in \fig{fig:single-gluon}, and therefore has no contribution to two-gluon correlator, i.e., ~\eq{eq:24}.

The calculations of the eight connected diagrams is similar to that of single gluon spectrum, except that  $f^{abc}f^{abd}=N_c\delta^{cd}$ and $f^{ade}f^{bef}f^{cfd}=\frac{N_c}{2}f^{abc}$. Four of the eight connected diagrams are shown in \fig{fig:3}, they give dominant contributions to the double-gluon inclusive production, and can be written in a compact form,
\begin{figure}[hbt]
\begin{centering}
\begin{tabular}{ccccccc}
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-1-2.eps}}}$&+&
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-1-3.eps}}}$&+&
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-2-1.eps}}}$&+&
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-3-1.eps}}}$\\
(a)& &(b)& &(c)& &(d)
\end{tabular}
\par\end{centering}
\caption{Four of the eight connected diagrams of double-gluon inclusive production in a collision. \fig{fig:3} (a) and (c) are called as the diffractive diagrams since the produced gluon are emitted from the same quark line in amplitude, \fig{fig:3} (b) and (d) are the interference diagrams.}
\label{fig:3}
\end{figure}

\begin{equation}\label{eq:18}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=\frac{\alpha^2_s N^2_c S_\perp}{\pi^{8}(N^2_c-1)^3}
\frac{1}{\bo{p}^2_{\perp}\bo{q}^2_{\perp}}
	\int \frac{d^2\bo{k}_\perp}{(2\pi)^2}(D_A+D_B),
\end{equation}
where
\begin{eqnarray}\label{eq:19}
D_A= \Phi^2_{A_1}(y_p,\bo{k}_\perp)\Phi_{A_2}(y_p,\bo{p}_\perp-\bo{k}_\perp)D_2,\nonumber\\
D_B= \Phi^2_{A_2}(y_q,\bo{k}_\perp)\Phi_{A_1}(y_p,\bo{p}_\perp-\bo{k}_\perp)D_1,
\end{eqnarray}
with
\begin{equation}
D_{1(2)}=\Phi_{A_{1(2)}}(y_q,\bo{q}_\perp+\bo{k}_\perp)
+\Phi_{A_{1(2)}}(y_q,\bo{q}_\perp-\bo{k}_\perp).
\end{equation}
It should be noticed that the pre-factor in \eq{eq:18} is consistent with that in ~\cite{Dusling2012b}, which is corrected from that in ~\cite{Dusling2009}.

\begin{figure}[hbt]
\begin{centering}
\begin{tabular}{ccccccc}
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-2-2.eps}}}$&+&
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-3-3.eps}}}$&+&
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-2-3.eps}}}$&+&
$\vcenter{\hbox{\includegraphics[scale=0.50]{figs/double-3-2.eps}}}$\\
(a)& &(b)& &(c)& &(d)
\end{tabular}
\par\end{centering}
\caption{Four of the eight connected diagrams for double-gluon inclusive production in a collision. \fig{fig:3} (a) and (c) are the double diffractive and non-diffractive diagrams, \fig{fig:3} (b) and (d) are the diffractive diagrams.}
\label{fig:4}
\end{figure}

The left four connected diagrams are shown in ~\fig{fig:4}, they give sub-leading contributions of $p_\perp/Q_s$ to double-gluon inclusive production, and can be written as,
\begin{equation}\label{eq:21}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=\frac{\alpha^2_s N^2_c S_\perp}{\pi^8 (N^2_c-1)^3}
\frac{1}{\bo{p}^2_\perp \bo{q}^2_\perp}
\sum_{j=\pm}
\left[D_3(\bo{p}_\perp,j\bo{q}_\perp)+
\frac{1}{2}D_4(\bo{p}_\perp,j\bo{q}_\perp)\right],
\end{equation}
where $D_3(\bo{p}_\perp,j\bo{q}_\perp)=\delta^2(\bo{p}_\perp - j\bo{q}_\perp)\left[\mathcal{I}^2_1 + \mathcal{I}^2_2 + 2\mathcal{I}^2_3\right]$, with
\begin{eqnarray}\label{eq:d3}
\mathcal{I}_1=
\int\frac{d^2\bo{k}_\perp}{(2\pi)^2}
\Phi_{A_1}(y_p,\bo{k}_{1\perp})\Phi_{A_2}(y_q,\bo{p}_\perp-\bo{k}_{1\perp})
\frac{(\bo{k}_{1\perp}\cdot\bo{p}_\perp - \bo{k}^2_{1\perp})^2}{\bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2},\nonumber\\
\mathcal{I}_2=
\int\frac{d^2\bo{k}_{1\perp}}{(2\pi)^2}
\Phi_{A_1}(y_p,\bo{k}_{1\perp})\Phi_{A_2}(y_q,\bo{p}_\perp-\bo{k}_{1\perp})
\frac{\vert\bo{k}_{1\perp}\times\bo{p}_\perp\vert^2}{\bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2},\nonumber\\
\mathcal{I}_3=
\int\frac{d^2\bo{k}_{1\perp}}{(2\pi)^2}
\Phi_{A_1}(y_p,\bo{k}_{1\perp})\Phi_{A_2}(y_q,\bo{p}_\perp-\bo{k}_{1\perp})
\frac{(\bo{k}_{1\perp}\cdot\bo{p}_\perp)\vert\bo{k}_{1\perp}\times\bo{p}_\perp\vert}{\bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2},
\end{eqnarray}

and
\begin{equation}\label{eq:d4}
\begin{split}
D_4(\bo{p}_\perp,j\bo{q}_\perp)=&
\int\frac{d^2\bo{k}_{1\perp}}{(2\pi)^2}
\Phi_{A_1}(y_p,\bo{k}_{1\perp})
\Phi_{A_1}(y_p,\bo{k}_{2\perp})
\Phi_{A_2}(y_q,\bo{p}_\perp-\bo{k}_{1\perp})
\Phi_{A_2}(y_q,\bo{p}_\perp-\bo{k}_{2\perp})\\
&\times
\frac{
(\bo{k}_{1\perp}\cdot\bo{p}_\perp-\bo{k}^2_{1\perp})
(\bo{k}_{2\perp}\cdot\bo{p}_\perp-\bo{k}^2_{2\perp})
+(\bo{k}_{1\perp}\times\bo{p}_\perp)\cdot(\bo{k}_{2\perp}\times\bo{p}_\perp)}{\bo{k}^2_{1\perp}(\bo{p}_\perp-\bo{k}_{1\perp})^2}\\
&\times
\frac{
(\bo{k}_{1\perp}\cdot j\bo{q}_\perp-\bo{k}^2_{1\perp})
(\bo{k}_{2\perp}\cdot j\bo{q}_\perp-\bo{k}^2_{2\perp})
+(\bo{k}_{1\perp}\times\bo{q}_\perp)\cdot
(\bo{k}_{2\perp}\times\bo{q}_\perp)}
{\bo{k}^2_{2\perp}(j\bo{q}_\perp-\bo{k}_{1\perp})^2},
\end{split}
\end{equation}
where $\bo{k}_{2\perp}=\bo{p}_\perp+j\bo{q}_\perp-\bo{k}_{1\perp}$.

\eq{eq:18} and \eq{eq:21} are corresponding to the results of \fig{fig:3} and \fig{fig:4}, respectively. Where \eq{eq:21} is suppressed by additional powers of $\bo{p}_\perp$ and $\bo{q}_\perp$ relative to the leading terms in \eq{eq:18}.

From the discussions above, the double-gluon inclusive production of 9 distinct diagrams can be identically written as,
\begin{equation}\label{eq:22}
\begin{split}
\left\langle\frac{dN_2}{d^2\bo{p}_\perp dy_p d^2\bo{q}_\perp dy_q}\right\rangle_{\mathrm{LLog}}
=\frac{\alpha^2_s N^2_c S_\perp}{\pi^{8}(N^2_c-1)^3}
\frac{1}{\bo{p}^2_{\perp}\bo{q}^2_{\perp}}
&	\left[\int \frac{d^2\bo{k}_\perp}{(2\pi)^2}(D_A+D_B)\right.\\
&\left.+\sum_{j=\pm}
\left[D_3(\bo{p}_\perp,j\bo{q}_\perp)+
\frac{1}{2}D_4(\bo{p}_\perp,j\bo{q}_\perp)\right]\right]\\
&+\left\langle\frac{dN_1}{d^2\bo{p}_\perp dy_p}\right\rangle_\mathrm{LLog}
\left\langle\frac{dN_1}{d^2\bo{q}_\perp dy_q}\right\rangle_\mathrm{LLog}.
\end{split}
\end{equation}
Where the first and second lines of ~\eq{eq:22} corresponds to ~\fig{fig:3} and \fig{fig:4}, respectively.  These are the main parts of two-gluon correlations. The third line of \eq{eq:22} corresponds to \fig{fig:1-1}, which has no contribution to two-gluon correlations.
\begin{figure}[hbt]
\begin{centering}
\begin{tabular}{ccc}
$\vcenter{\hbox{\includegraphics[scale=0.60]{figs/bfkl.eps}}}$&\hspace{1cm} &
$\vcenter{\hbox{\includegraphics[scale=0.60]{figs/short-range.eps}}}$\\
(a)& &(b)
\end{tabular}
\par\end{centering}
\caption{The amplitude of two-gluon production in a collision, the left diagram (a) gives back-to-back correlations, and the right diagram (b) gives near-side collinear correlations.}
\label{fig:other}
\end{figure}

In addition, according to the standard QCD approach, the double gluon inclusive production should still contain other two Feynman diagrams as showed in \fig{fig:other}. Where \fig{fig:other}(a) is the dynamics of Mueller-Navelet jets~\cite{Fadin1997, Colferai2010}. It gives the dominant contribution to away-side $\Delta\phi=\pi$ correlations and is negligible in the near-side $\Delta\phi=0$ correlations. \fig{fig:other}(b) is gluon jets stem from collinear jet shower. So they are both out of the concerning region of CGC dynamics.

\section{Two-gluon rapidity correlations in the leading and sub-leading orders}
\label{sec:results}

The values of the sub-leading order (\eq{eq:21}) and how it is relative to that of the leading order (\eq{eq:18}) are concerned and interested. In this section, we will discuss their numerical results.

\begin{figure}[t]
\includegraphics[scale=0.6]{figs/trans.pdf}
\caption{(a)Transverse momentum dependence of two-gluon rapidity correlations in fixed rapidity gap in $pp$ collisions at $\sqrt{s}=7\ \mathrm{TeV}$.  Where $y_p=0,\Delta\phi=0$. The leading and sub-leading diagrams are presented by the solid and dashed lines, respectively.
Three rapidity gaps are $\Delta y=0,1.5$ and 3, and presented by the red, blue and black lines, respectively. (b) The ratios of sub-leading order to that of leading ones for the same rapidity gaps as (a).}
\label{pp-trans}
\end{figure}

To find the transverse momentum region for the strongest correlations as showed in our previous paper~\cite{Zhao2015}, we present firstly in ~\fig{pp-trans}(a) the transverse momentum dependence of two-gluon correlation in given rapidity gaps in $pp$ collisions at $\sqrt{s}=7\ \mathrm{TeV}$ with running coupling.  Where the relative azimuthal angle between two gluon is set to zero i.e., $\Delta\phi=0$, as the correlation is the maximum when the transverse momentum of two selected gluon is collimation, i.e., $\Delta\phi=\phi_q-\phi_p=0,\pi$~\cite{Dumitru2011d,Dusling2012b}. For simplicity, the transverse momentum of two gluon is set to $p_\perp=q_\perp$ and $y_p=0$ is fixed. The leading and sub-leading contributions are presented by the solid and dashed lines, respectively. Three fixed rapidity gaps are $\Delta y\equiv y_{q}-y_{p}=0,1.5$ and 3 presented by the red, blue and black lines, respectively.

\fig{pp-trans} (a) shows again that the transverse momentum dependence of the two-gluon rapidity correlations is peaked around the summation of saturation momenta of two incident particles, i.e., $p_\perp\sim Q_{sA}+Q_{sB}$,  consistent with our previous results with fixed coupling constant~\cite{Zhao2015}. However, the correlation with running coupling is much stronger than that with fixed coupling constant, cf., the ~Fig.3 of~\cite{Zhao2015}.

In contrast, the contributions of sub-leading order are much smaller than that of leading ones, and almost independent of the transverse momentum. The latter is because of the contributions of sub-leading order in ~\eq{eq:d3} and ~\eq{eq:d4} is further integrated to the transverse momentum of the original uGD, and therefore smeared the influence of saturation momentum. This shows that the contributions of sub-leading order are not as sensitive as that of the leading ones to the effect of gluon saturation.

From the transverse momentum dependence of two-gluon rapidity correlations at three different rapidity gaps, it is clear that the correlation decreases with the width of the rapidity gap. The smallest rapidity gap, the red dash line, has the strongest correlation. The contributions of the sub-leading order are mainly short-range rapidity correlations.

To see the relative contribution of sub-leading order, we present in the~\fig{pp-trans} (b), the ratios of sub-leading order to that of the summation of the leading and sub-leading ones in three same rapidity gaps as in ~\fig{pp-trans} (a). Where the red, blue and black dashed lines correspond to $\Delta y=0,1.5$ and 3, respectively.

The ~\fig{pp-trans} (b) shows that the contributions of sub-leading order is less than $10\%$ for larger rapidity gap, such as $\Delta y=3$, showed by solid black line in~\fig{pp-trans} (b), and larger than $10\%$ for small rapidity gap, such as $\Delta y\le 1.5$ or smaller, showed by solid blue and red lines in~\fig{pp-trans} (b). So the contribution of sub-leading order is important in short-range rapidity correlations.

\begin{figure}[hbt]
\includegraphics[width=\linewidth]{figs/rapidity.pdf}
\caption{Two-gluon rapidity correlations in the leading (a) and sub-leading (b) orders for $pp$ collisions at $\sqrt{s}=7$ TeV for three sets of transverse momenta ($\bo{p}_\perp=\bo{q}_\perp=1$ (red) ,2(blue) and $3$ GeV (black), respectively. Where (c) is the summation of (a) and (b), and (d) is the ratio of (b) to (c). }
\label{fig:pp-rap}
\end{figure}

Now we turn to study the two-gluon rapidity correlations around transverse momentum $p_\perp\sim Q_{sA}+Q_{sB}$. For fixed transverse momentum set $\bo{p}_\perp=\bo{q}_\perp$, the two-gluon rapidity correlations in the leading and sub-leading orders for $pp$ collisions at $\sqrt{s}=7$ TeV are presented in ~\fig{fig:pp-rap}(a) and 7(b), respectively. Where again, $y_p=0$, $\Delta\phi=0$. Three sets of transverse momenta are $\bo{p}_\perp=\bo{q}_\perp=1$ (red line) , 2(blue line) and $3$ GeV (black line), respectively.

\fig{fig:pp-rap}(a) shows again that the strongest correlation happens when the transverse momenta of two gluon are both the summation of saturation momenta of two incident particles, i.e., $\bo{p}_\perp=\bo{q}_\perp\sim Q_{sA}+Q_{sB}=2$GeV, as showed by blue line. For each set of transverse momentum, the correlation is enhanced in both short and long rapidity ranges, in contrast to our previous results with fixed coupling constant~\cite{Zhao2015}, where the correlation is only enhanced in long rapidity range, and short-range correlations are almost flat, cf., Fig.~ 4 of ~\cite{Zhao2015}.

This is understandable. In the framework of CGC EFT, short-range rapidity correlation comes from the quantum correlation. When the coupling strength decreases with transverse momentum, i.e., the running coupling, the short-range correlation caused by the quantum correlations will change with coupling strength as well. The smaller the transverse momentum, the larger the coupling and the stronger correlation are. Therefore, the short-range rapidity correlation between two gluon with smaller transverse momentum is enhanced in the case of running coupling.

\fig{fig:pp-rap}(b) shows that the contributions of the sub-leading order are mainly short-range rapidity correlations, and the smaller the transverse momentum set has the stronger correlations. So short-range rapidity correlation caused by quantum correlation keeps in the sub-leading order as well as that of the leading one. Where the ridge-like long-range rapidity correlation due to the gluon of colour source is negligible in the sub-leading order.

To see the total contributions of the leading and sub-leading orders, we present in \fig{fig:pp-rap} (c) the summation of the leading order, the \fig{fig:pp-rap}(a),  and the sub-leading order, the
\fig{fig:pp-rap}(b). It shows again that the sub-leading order gives visible contributions to short-range rapidity correlations. To quantify the relative contribution of sub-leading order, in \fig{fig:pp-rap}(d), we present the ratio of sub-leading order to the summation one, i.e., \fig{fig:pp-rap}(c). It shows that when rapidity gap $|\Delta y|\le 2$, the contribution of the sub-leading order is larger than $10\%$. So the contribution of sub-leading order is not negligible in short-range rapidity correlation.

For $pA$ and $AA$ collisions,  the similar results are obtained. Where only the ratio of $AA$ is  slightly higher. This is because the saturation scale of the nucleus is larger than that of a proton at the same incident energy. It makes the contributions of both leading and sub-leading orders smaller, and so a slightly higher ratio.

So in the framework of CGC EFT, the stronger correlations in both short and long rapidity ranges are naturally obtained in $pp$, $pA$ and $AA$ collisions. This is qualitatively consistent with experimental data of $pPb$ and $PbPb$ collisions published by LHC collaboration~\cite{CMSCollaboration2012a,Chatrchyan2011},  where when the transverse momentum of the triggered and associated particles are in the middle region, i.e., $1\le p_{\perp}\le 3$ GeV/c, and relative azimuthal angle is near-side, i.e., $\Delta \phi=0$, the stronger two-particle correlations in short- and long-pseudo-rapidity ranges are observed.

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

In the framework of Color Glass Condensate (CGC) effective field theory (EFT), we calculate two-gluon rapidity correlations in the leading and sub-leading orders with running coupling constant.

The transverse momentum dependence of two-gluon rapidity correlation at three rapidity gaps are firstly presented. It is shown again that in the leading order, the correlation is still peaked at the summation of two saturation momenta of incident particles, the same as our previous results with fixed coupling~\cite{Zhao2015}. In contrast, the contributions of sub-leading order are almost independent of transverse momentum. So the contribution of the leading order is sensitive to the saturation momentum of incident particles, but the sub-leading one is not. Meanwhile, the results of different rapidity gaps indicate that the short-range rapidity correlation is important in the sub-leading order.

Then, the rapidity correlation pattern in the leading order show that both short and long-range rapidity correlations are enhanced, in contrast to our previous results with fixed coupling constant, where only ridge-like long-range rapidity correlations are obtained. Therefore, the running coupling enhances the short-range rapidity correlations. While, the contribution of sub-leading order is much smaller than that of leading one, and mainly short-range rapidity correlations.

As we pointed out in our previous work~\cite{Zhao2015}, the ridge-like long-range rapidity correlations is due to the colour source gluon. It is essential in the leading order, but negligible in sub-leading one. While, the short-range rapidity correlations come from quantum correlation. This quantum effect is same important in both leading and sub-leading orders.

Although the contributions of sub-leading order is much smaller than that of leading one,
its contributions to the short-range rapidity ($\Delta y\le 2$) correlation are larger than $10\%$, and therefore not negligible.

In the framework of CGC EFT, the stronger short- and long-range rapidity correlations are naturally obtained. This is qualitatively consistent with LHC experimental data on two-particle correlations~\cite{CMSCollaboration2012a,Chatrchyan2011}. So CGC EFT provides a good understanding for the features of two-particle rapidity correlations in high energy collisions.

\section*{Acknowledgments}
This work is supported in part by the Major State Basic Research Development Program of China under Grant No. 2014CB845402, the NSFC of China under Grant No. 11521064,
and the Ministry of Scicence and Technology (MoST) under grant No. 2016YFE0104800.