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

    $W^{\pm}$ vector bosons are heavy particles with masses of $ m_{W^\pm} 
\approx 80$~GeV/$c^{2}$~\cite{ParticleDataGroup:2020ssz,CDF:2022hxs}. They are 
mainly produced in the hard partonic scattering processes with large momentum 
transfer at the early stage of the (ultra-)relativistic nuclear 
collisions. In comparison with evolution time of the heavy-ion collision 
system, 10 to 100~fm/$c$ for instance, the $W^{\pm}$ decay time 
($\sim 0.0922~{\rm fm/}c$, estimated with full decay width, $t=h/\Gamma$ 
\cite{ParticleDataGroup:2020ssz}) is very short. The $W^{\pm}$ leptonic decays 
\begin{equation*}
    W^{\pm}\rightarrow l^{\pm}\nu_{l},~(l{\rm :}~e,\mu,\tau)
\end{equation*}
are nearly instantaneous. As the produced leptons weakly interact with the 
partonic and hadronic matter created in nuclear collisions, the $W^{\pm}$ is a 
powerful probe for investigating the properties of the initial state in the 
nuclear collision system. The main production processes of $W^\pm$ are 
\begin{equation*}
    u\bar{d} \rightarrow W^{+}, \hspace{0.4cm} d\bar{u} \rightarrow W^{-}
\end{equation*}
in the leading-order (LO) approximation~\cite{ParticleDataGroup:2020ssz,Martin:1999ww}.
With the involvement of the valence $u$- and $d$-quarks inside the nuclear 
system, the isospin difference between proton and neutron is important for the 
production of $W^{+}$ and $W^{-}$. It is therefore expected that the $W^\pm$ 
boson charge asymmetry and the corresponding decayed lepton charge asymmetry 
are related to the relative abundance of protons and neutrons (to the relative 
abundance of valence $u$- and $d$-quarks) of the collision system. The lepton 
charge asymmetry $A_{l}$ could be defined as the difference between the 
multiplicity of positive and negative leptons divided by their sum:
\begin{equation}
    A_{l} = \frac{N_{l^{+} \leftarrow W^{+}} - N_{l^{-} \leftarrow W^{-}}}
	    {N_{l^{+} \leftarrow W^{+}} + N_{l^{-} \leftarrow W^{-}}}.
\label{eq:A_l}
\end{equation}
Similarly, the collision system valence $u$ ($d$) quark number asymmetry 
$A_v$ is defined as
\begin{equation}
A_v = \frac{N_u - N_d}{N_u + N_d},
\end{equation}
where $N_u$ ($N_d$) refers to the number of valence quark $u$ ($d$) in 
the nuclear collision system.

    The lepton charge asymmetry have been measured in $p\bar p$ collisions at 
$\sqrt{s}=1.8$ and $1.96$~TeV by CDF~\cite{CDF:1991igh,CDF:1994cas,CDF:1998uzn,CDF:2005cgc,CDF:2009cjw} 
and D0~\cite{D0:2007pcy,D0:2008cgv,D0:2013xqc} collaborations on Fermilab 
Tevatron. They give early complementary constrains on Parton Distribution 
Function (PDF) relative to the deep inelastic scattering experiments. 
More precise measurements in $pp$ collisions at $\sqrt{s}$=7 and 8~TeV are 
reported by ATLAS~\cite{ATLAS:2011pph,ATLAS:2011qdp,ATLAS:2016nqi,ATLAS:2019fgb}, 
CMS~\cite{CMS:2011aa,CMS:2012ivw,CMS:2012fgk,CMS:2013pzl,CMS:2016qqr} and 
LHCb~\cite{LHCb:2012lki,LHCb:2014liz,LHCb:2015mad,LHCb:2016zpq} collaborations 
on the Large Hadron Collider (LHC) at CERN. ATLAS also give the measurements 
in $pp$ collisions at $\sqrt{s}=2.76$~\cite{ATLAS:2019fyu} and $5.02$~TeV
~\cite{ATLAS:2018pyl}. Beyond the elementary $pp$ collisions, the lepton charge 
asymmetry has been measured in the p+Pb collisions at 
$\sqrt{s_{\rm NN}}=5.02$~TeV by ALICE~\cite{ALICE:2016rzo} and at 
$\sqrt{s_{\rm NN}}=8.16$~TeV by ALICE~\cite{ALICE:2022cxs} and CMS~\cite{CMS:2019leu}. 
The results of lepton charge asymmetry in Pb+Pb collisions at 
$\sqrt{s_{\rm NN}}=2.76$~TeV are first given by ATLAS~\cite{ATLAS:2014sic} and 
CMS~\cite{CMS:2012fgk}. Recently, ATLAS and ALICE published their results 
of $W^\pm$ production cross-section and lepton charge asymmetry in Pb+Pb 
collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV~\cite{ATLAS:2019ibd,ALICE:2022cxs}. 
All those measurements are declared to be well reproduced by the LO and/or 
next-to-leading-order (NLO) perturbative Quantum ChromoDynamics (pQCD) 
calculations~\cite{Dulat:2015mca,Eskola:2009uj,Kusina:2016fxy,Eskola:2016oht}
using the CT14 PDF set~\cite{Dulat:2015mca} with or without the nuclear 
modified PDF (nPDF).

    In this paper, the parton and hadron cascade model PACIAE \cite{Sa:2011ye} is 
employed to reproduce the ATLAS and ALICE lepton charge asymmetry data 
in Pb+Pb at $\sqrt{s_{\rm NN}}=5.02$~TeV~\cite{ATLAS:2019ibd,ALICE:2022cxs}. We 
then calculate the lepton charge asymmetry as a function of the collision 
system valence quark number asymmetry in both the elementary nucleon-nucleon 
(NN) and the nucleus-nucleus (AA) collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV. An 
interesting linear scaling behavior in the lepton charge asymmetry as the 
function of collision system valence quark number asymmetry is observed 
among the different NN and AA collision systems.

    The paper is organized as follows. In Sec.~\ref{sec:mod}, we describe the 
PACIAE model and the setups for $W^\pm$ production in \textsc{Pythia} 6.4, 
with a precision equivalent to the NLO pQCD. The simulated results and 
discussions are presented in Sec.~\ref{sec:res}. At last, a summary is 
given in Sec.~\ref{sec:summ}.

\section{\label{sec:mod} Model}

    A parton and hadron cascade model PACIAE~\cite{Sa:2011ye} is employed to 
simulate the elementary nucleon-nucleon (NN: $pp$, $pn$, $np$ and $nn$) 
and the nucleus-nucleus (AA: Cu+Cu, Xe+Xe, Au+Au, Pb+Pb and U+U) 
collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV.

    The PACIAE model is based on \textsc{Pythia} event generator (version $6.4.28$) 
\cite{Sjostrand:2006za}. For NN collisions, with respect to \textsc{Pythia}, 
the partonic and hadronic rescatterings are introduced before and after 
the hadronization, respectively. The final hadronic state is developed 
from the initial partonic hard scattering and parton showers, followed by parton
rescattering, string fragmentation, and hadron rescattering stages. Thus, 
the PACIAE model provides a multi-stage transport description on the evolution 
of the NN collision system.

    For AA collisions, the initial positions of nucleons in the colliding nucleus 
are sampled according to the Woods-Saxon distribution. Together with the 
initial momentum setup of $p_{x} = p_{y} = 0$ and $p_{z} =p_{\rm beam}$ for 
each nucleon, a list containing the initial state of all nucleons in a given 
AA collision is constructed. A collision happened between two nucleons 
from different nuclei if their relative transverse distance is less than or 
equal to the minimum approaching distance:
$D\leq\sqrt{\sigma_{\rm NN}^{\rm tot}/\pi}$. The collision time is calculated 
with the assumption of straight-line trajectories. All such nucleon pairs 
compose an NN collision time list. An NN collision with least collision time is 
selected from the list and executed by \textsc{Pythia} (PYEVNW subroutine) with 
the hadronization temporarily turned-off, as well as the strings and diquarks 
broken-up. The nucleon list and NN collision time list are then updated. A new 
NN collision with least collision time is selected from the updated NN collision 
time list and executed by \textsc{Pythia}. With repeating the aforementioned 
steps till the NN collision list empty, the initial partonic state is 
constructed for a AA collision.

    Then, the partonic rescatterings are performed, where the LO-pQCD parton-parton 
cross section~\cite{Combridge:1977dm,field} is employed. After partonic 
rescattering, the string is recovered and then hadronized with the Lund string 
fragmentation scheme resulting in an intermediate hadronic state. Finally, the 
system proceeds into the hadronic rescattering stage and produces the final 
hadronic state observed in the experiments.

    The PACIAE Monte-Carlo simulation provides a complete description of the NN 
and AA collisions, which includes the partonic initialization, partonic 
rescattering, hadronization, and the hadronic rescattering stages. Meanwhile, 
a PACIAE model simulation could be selected to stop at any stage desired 
conveniently. In this work, the simulations are stopped at the partonic 
initialization stage, after the initial hard scattering and parton shower 
processes happened.

    In this study, the $W^\pm$ relevant production channels are activated in a 
user-controlled approach by setting ${\rm MSEL}=0$ in $\textsc{Pythia}$, 
together with the following subprocesses switched-on:
\begin{align*}
    f_{i}\overline{f}_{j} & \rightarrow W^{+}/W^{-},   \\
    f_{i}\overline{f}_{j} & \rightarrow g W^{+}/W^{-}, \\
    f_{i}\overline{f}_{j} & \rightarrow \gamma W^{+}/W^{-}.
\end{align*}
In the above equations, $f$ refers to fermions (quarks) and its subscript 
stands for the flavor code.

    In addition, the parton shower is switched off by setting MSTP(61)=0 and 
MSTP(71)=0, to prevent the double-counting. 
The PDF and other parameters are kept the same as the default \textsc{Pythia}. 
Thus the PACIAE (\textsc{Pythia}) model simulation is corresponding to the NLO 
pQCD in precision.

\section{\label{sec:res} Results and Discussions}

    In the left panel of Fig.~\ref{fig:comp}, we compare the PACIAE model 
calculated electron and muon combined lepton charge asymmetry (black solid 
squares) with ATLAS data~\cite{ATLAS:2019ibd} (red solid circles) for the 
0-80\% centrality class in Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV. 
The bottom sub-panel shows the differences between PACIAE calculation and ATLAS 
data, where the shaded boxes show the total uncertainties in experimental data. 
The kinematic cuts in PACIAE simulation is the same as the ones in ATLAS 
experiments~\cite{ATLAS:2019ibd}, i.e.
\begin{align*}
    p^{l}_T > 25 {\rm~GeV}/c, & \hspace{0.4cm} p^{\nu}_T > 25 {\rm~GeV}/c, \\
    |\eta_l| < 2.5, & \hspace{0.4cm} m_T > 40 {\rm~GeV}/c^2,
\end{align*}
where $p^{l}_T$ and $p^{\nu}_T$ stand for the transverse momentum of the 
lepton and the corresponding neutrino, respectively. 

\begin{figure*}[htbp]
    \includegraphics[width=0.8\textwidth]{PACIAE_ATLAS_ALICE.eps}
    \caption{\label{fig:comp} Left panel: the combined electron and muon 
charge asymmetry as a function of absolute pseudorapidity for the 0-80\% 
centrality class in Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV. The black 
solid squares with vertical error bar (hard to see) show the results of PACIAE 
simulations. The red solid circles with vertical error bars and shaded error 
boxes represent the ATLAS data ~\cite{ATLAS:2019ibd}. The lower sub-panel shows 
the differences between the PACIAE calculations and the ATLAS data, with the 
shaded boxes representing the total experimental uncertainties. The 
PACIAE result points are shifted horizontally for better visibility.
Right panel: the muon charge asymmetry for the 0-90\% centrality class in 
Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV, with kinematic cuts of 
$p^{\mu}_T > 10$~GeV/$c$ and $2.5 < y_{\mu} < 4.0$. Here the open square with 
vertical error bar shows the result of PACIAE simulation, and the blue solid 
circle with vertical error bar and shaded error box represents ALICE 
darum~\cite{ALICE:2022cxs}.}
\end{figure*}

\noindent $|\eta_l|$ refers to the
absolute pseudorapidity of the lepton, and $m_T$ is the transverse mass of the 
lepton and neutrino system. 
The pseudorapidity bins are divided as follows:
\begin{align*}
    |\eta_l| = &~(0,0.21),~~~~~(0.21,0.42), ~(0.42,0.63), ~(0.63,0.84), \\
               &~(0.84,1.05), ~(1.05,1.37), ~(1.37,1.52), ~(1.52,1.74), \\
               &~(1.74,1.95), ~(1.95,2.18), ~(2.18,2.5).
\end{align*}

    We see in this panel, 
the PACIAE calculated lepton charge asymmetry is well consistent with ATLAS 
data~\cite{ATLAS:2019ibd}, except the most forward $|\eta_l|$ bin. This 
exception is also existed in the pQCD predictions calculated with the CT14 NLO 
PDF set, even in the EPPS16 and nCTEQ15 nPDF sets~\cite{ATLAS:2019ibd}.

    The PACIAE calculation of the muon charge asymmetry (open square) for the 
0-90\%  centrality class in Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV is 
compared with the ALICE datum ~\cite{ALICE:2022cxs} (blue solid circle) in the 
right panel of Fig.~\ref{fig:comp}. The same kinematic cuts of 
$p^{\mu}_T > 10$~GeV/$c$ and $2.5 < y_{\mu} < 4$~\footnote[1]
{The ALICE muon spectrometer covers negative pseudorapidity. However, in 
symmetric Pb+Pb collisions, positive values of rapidity are conventionally
used for the muon coverage.} 
are set in both the experiment and theory. One can see in this panel that, the 
calculated result agrees well with the ALICE datum. 

    In Fig.~\ref{fig:Amu_pp}, we give the muon charge asymmetry as a function of 
the valence quark number asymmetry in the elementary NN collisions of $pp$ 
(open triangle-ups), $pn$ (triangle-lefts), $np$ (triangle-rights), and $nn$ 
(triangle-downs) 

\begin{figure}[htbp]
    \includegraphics[width=0.45\textwidth]{Amu_pp_vs_Av.eps}
    \caption{\label{fig:Amu_pp} The muon charge asymmetry as a function of the 
valence quark asymmetry in the elementary NN collisions ($pp$, open 
triangle-ups; $pn$, open triangle-lefts; $np$, open triangle-rights; 
and $nn$, open triangle-downs) at $\sqrt{s}=5.02$~TeV for three rapidity 
intervals of $|y|<1$, $2.5<|y|<4.0$, and full rapidity phase-space. The dashed, 
dotted, and dash-dotted lines are the corresponding linear fitting to the $A_\mu$ 
vs. $A_v$ functions, respectively. The kinematic cut is $p^\mu_T > 10$~GeV/$c$. 
The red short-dashed zero-baseline and the shaded area are added for better 
distinguishment between the positive and negative values of $A_\mu$.}
\end{figure}

\noindent at $\sqrt{s}=5.02$~TeV. The kinematic cut of $p^\mu_T > 
10$~GeV/$c$ is implemented in the simulations for three $y_\mu$ intervals: 
$|y_\mu|<1$, $2.5<|y_\mu|<4.0$, and full $y_\mu$ phase-space. 
The dashed, dotted and dash-dotted lines represent the linear fitting to the 
function of $A_\mu$ vs. $A_v$ in three $y_\mu$ intervals above, respectively. 
One can see all of $A_\mu$ in $pp$, $pn$, $np$ and $nn$ collisions are positive 
in mid-rapidity ($|y_\mu| < 1$), while negative in the forward rapidity 
($2.5 < |y_\mu| <4$). This sign-changing phenomenon has been also observed in 
LHCb measurements~\cite{LHCb:2012lki,LHCb:2014liz,LHCb:2015mad,LHCb:2016zpq}. 
It comes from the different helicity dependence of the lepton couplings to the 
boson~\cite{LHCb:2012lki}. For full rapidity phase-space, $A_\mu$ in $pn$ and 
$np$ collisions is consistent with zero owing to the vanishing relative numbers 
between $u$- and $d$-valence quarks, while the $A_\mu$ in $pp$ and $nn$ are 
almost equal in magnitude but opposite in sign. On the other hand, we can see 
the slope fitted in mid-rapidity is gentler than the one in the forward 
rapidity, due to the isospin effect in $W^{\pm}$ production is more pronounced 
at forward rapidity~\cite{ConesaDelValle:2007dza}.


    The correlation between $A_\mu$ and $A_v$ in symmetrical AA collisions of Cu+Cu 
(solid circles), Xe+Xe (solid squares), Au+Au (solid stars), Pb+Pb (solid 
pentagons) and U+U (solid diamonds) at $\sqrt{s_{\rm NN}}=5.02$~TeV for three 
rapidity intervals is shown in Fig.~\ref{fig:Amu_AA}, with the kinematic cut of 
$p^\mu_T > 10$~GeV/$c$. The corresponding linear fitting is also given. 
In this figure, the valence quark number asymmetry ($A_v$) are all negative, 
since the number of neutrons (the number of $d$ quarks) is always greater than 
that of protons ($u$ quarks) in heavy nuclei. The $A_\mu$ is positive 
in the mid-rapidity bin for all five AA collision systems but negative 
in the forward and full rapidity regions. It can be attributed to the $V-A$ 
decays of $W$ boson. This structure also affects the slope of the linear 
fitting. With the ratio of neutron to proton ($A_v$ magnitude) increasing, the 
fraction of $W^- \rightarrow l^- \nu_l$ becomes larger in forward rapidity 
resulting a steeper slope than the one in mid-rapidity~\cite{LHCb:2012lki,ConesaDelValle:2007dza}. 
It is more clear in Fig.~\ref{fig:Amu_AApp}, where we draw the $A_\mu$ vs. $A_v$ 
function in all the NN and AA collisions together. Here we see the fitted slope 
in the full rapidity is larger than the one in mid-rapidity, but is smaller than 
that in forward rapidity. One more interesting feature is that, the linear 
scaling behavior of $A_\mu$ vs. $A_v$ shown in Fig.~\ref{fig:Amu_pp} is also 
observed in Fig.~\ref{fig:Amu_AA} and~\ref{fig:Amu_AApp}.

\begin{figure}[htbp]
    \includegraphics[width=0.45\textwidth]{Amu_AA_vs_Av.eps}
    \caption{\label{fig:Amu_AA} The muon charge asymmetry as a function of 
the valence quark number asymmetry in the symmetric AA collisions (Cu+Cu, solid 
circles; Xe+Xe, solid squares; Au+Au, solid stars; Pb+Pb, solid pentagons; and 
U+U, solid diamonds) at $\sqrt{s_{\rm NN}}=5.02$~TeV for three rapidity 
intervals of $|y_\mu|<1$, $2.5<|y_\mu|<4.0$, and full rapidity phase-space. The 
dashed, dotted and dash-dotted lines are the corresponding linear fitting to the 
$A_\mu$ vs. $A_v$ functions, respectively. The kinematic cut is 
$p^\mu_T > 10$~GeV/$c$. The red short-dashed zero-baseline and the shaded area 
are added for better distinguishment between the positive and negative values 
of $A_\mu$.}

\end{figure}

\begin{figure}[htbp]
    \includegraphics[width=0.45\textwidth]{Amu_AApp_vs_Av.eps}
    \caption{\label{fig:Amu_AApp} The muon charge asymmetry as a function of 
the valence quark number asymmetry in the elementary NN collisions and the symmetric 
AA collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV, obtained by putting 
Fig.~\ref{fig:Amu_pp} and \ref{fig:Amu_AA} together.}
\end{figure}

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

    In summary, we studied the lepton charge asymmetry from the $W^\pm$ decays in 
both the NN and AA nuclear collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV in frame 
of the Monte-Carlo simulation model PACIAE. The simulated results well 
reproduce the corresponding ATLAS and ALICE data. Furthermore, the correlation 
between the lepton charge asymmetry $A_\mu$ and the valence quark number 
asymmetry $A_v$ is studied in the above nuclear collisions for three rapidity 
intervals: $|y_\mu|<1$, $2.5 < |y_\mu| < 4$, and the full rapidity phase space. 
In each rapidity interval, a linear scaling behavior of $A_\mu$ vs. $A_v$ 
function is observed in the elementary NN collision systems 
($pp$, $pn$, $np$, and $nn$), the AA collision systems 
(Cu+Cu, Xe+Xe, Au+Au, Pb+Pb, and U+U), and even in their combination.

    The universal linear scaling behavior between $A_\mu$ and $A_v$ explored in 
both nucleon-nucleon and nucleus-nucleus collision systems may provide an 
additional constraint in the extraction of PDF (nPDF). On the contrary, we 
also note that the nuclear modifications on the vacuum nucleon PDF change the 
in-going parton kinematics in the hard scattering processes and consequently 
impact the $p_T$- and $y$-distribution of the out-going particles. This effect 
should result in a breaking of the $A_\mu$ vs. $A_v$ scaling established if 
looking into a given lepton $p_T$--$y$ phase space. Therefore, by measuring 
$p_T$ and $y$ differentially, the deviation of the $A_\mu$ vs. $A_v$ 
correlations between NN and AA collision systems should provide a new tool to 
constrain the nPDF in a wide Bjorken-$x$ region. This has to be studied 
further.

\begin{acknowledgments}
    The authors thank helpful discussion with Ming-Rui Zhao. This work was 
supported by the National Natural Science Foundation of China (11775094, 
11805079, 11905188, 11775313, 11905163), the Continuous Basic Scientific 
Research Project (No.WDJC-2019-16) in CIAE, National Key Research and 
Development Project (2018YFE0104800) and by the 111 project of the foreign 
expert bureau of China.
\end{acknowledgments}

