Maximally entangled proton and charged hadron multiplicity in Deep Inelastic Scattering

We study the proposal by Kharzeev-Levin to determine entanglement entropy in Deep Inelastic Scattering (DIS) from parton distribution functions (PDFs) and to relate the former to the entropy of final state hadrons. We find several uncertainties in the current comparison to data, in particular the overall normalization, the relation between charged versus total hadron multiplicity in the comparison to experimental results as well as different methods to determine the number of partons in Deep Inelastic Scattering. We further provide a comparison to data based on leading order HERA PDF as well as PDFs obtained from an unintegrated gluon distribution subject to next-to-leading order Balitsky-Fadin-Kuraev-Lipatov and Balitsky-Kovchegov evolution. Within uncertainties we find good agreement with H1 data. We provide also predictions for entropy at lower photon virtualities, where non-linear QCD dynamics is expected to become relevant.


Introduction
Entanglement is a nonlocal correlation unique to quantum systems [1], see also the reviews [2,3]. There are various proposals how to study entanglement in high energy physics such as neutrino oscillations, spin correlations of t −t quarks or Λ hyperons [4][5][6][7]. A measure of entanglement which is of particular interest is entanglement entropy [8]. In [9] it has been proposed that one can associate entanglement entropy with the system of partons probed in Deep Inelastic Scattering (DIS) experiments. The proposal necessarily requires a measurement process which introduces a bi-partition to the system. The measured system is no longer in a pure quantum state and as a consequence entropy arises. The proposal in [9] considers coordinate space entropy and is motivated by exact results obtained in conformal field theories in 1+1 dimensions [10,11]. The 1+1 dimensional picture provides reasonable guidance, since the basic quantity which describes the system of partons at some resolution scale is the integrated density of partons, i.e. the parton distribution function (PDF). It provides information 1/Q 1/x S(x,Q) Figure 1: The figure illustrates the generation of entanglement in DIS within the dipole picture with Q the photon virtuality. The blue region represents the proton, the red circle the area singled out by the virtual photon. The segments represent color dipoles in color singlet states. The entanglement arises due to dipoles that are partially in the red and blue region. about the one dimensional spatial structure of protons, although the underlying dynamics may take place in transverse dimension as well. For further developments and other approaches see [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30].
The measurement process is provided by a probe. In DIS this probe is given by the virtual photon that is exchanged between electron and proton. Since the photon probes only parts of the partonic system of the proton, the remaining part has to be integrated out or traced over, giving rise to a reduced density matrix and therefore entanglement entropy. We note that this approach was also used to estimate entanglement entropy produced in proton-proton collisions at the LHC [13]. It can be motivated as follows: for DIS at low x and referring to the proton rest frame, it is possible to take color dipoles as the basic degrees of freedom. Color dipoles are color singlets and therefore natural candidates to generate entanglement entropy 1 . The mechanism of the entanglement is the following: • once the virtual photon resolves substructure in the proton, it singles out a certain region of area 1 See also [22] where the dipole degrees of freedom are used to obtain entanglement entropy ∼ 1/Q 2 see Fig. 1.
• it might happen that one of the quarks that is a constituent of the dipole is inside the singled out region A and the other part is in the region B. Since the dipole is a color singlet state, the quark-antiquark system is strongly correlated and provides the source of entanglement between both regions. Due to the interaction, the dipole breaks up and the quark in the region B gives rise to the observed final state hadrons.
• the other quark gives a contribution to the parton distribution function measured in the DIS process • since both constituents of the dipole were in different regions, the entropy of both regions is identical. Only dipoles that bridge the regions give rise to the observed entropy.
• as one goes to higher energies the number of dipoles increases (higher twists) and more and more dipoles are partially in region A and region B. As a consequence, entropy grows.
• turning to smaller Q 2 , more and more dipoles are completely inside region A and the entropy should eventually decrease.
The framework presented in [9] provides furthermore an explicit formula on how to calculate this entanglement entropy as well as means how to obtain the latter from data through the determination of the entropy of the observed final state hadrons. Since the entire approach is based on distributions of color dipoles, which themselves can be derived from QCD within high energy factorization with x 1 the expansion parameter, the proposed description is naturally restricted to the low x regime. NLO fits of low x gluon distribution are found to provide a good description of HERA data in the region x < 0.01, see [31,32] for a fit with the unintegrated gluon evolved with next-to-leading (NLO) Balitsky-Fadin-Kuarev-Lipatov (BFKL) evolution. The same phase space restriction x < .01 has been used in the fits based on leading order running coupling [33] and NLO [34] Balitsky-Kovchegov evolution, which is the non-linear low x evolution equation formulated in terms of the dipole amplitude. It is therefore natural to use x < 0.01 as an upper limit on phase space region, where the above arguments can be expected to be applicable. It is needless to say that more conservative bounds would yield even stronger justification for the underlying low x approximations.
In the recent paper [35] we provided numerical evidence for this proposal, see also [36] for related work. In this way entropy at the partonic level is determined through where g(x, µ 2 f ) denotes the gluon distribution function at the factorization scale µ f and Σ(x, µ 2 f ) = ∑ n f a=1 q a (x, µ 2 ) +q a (x, µ 2 ) the quark flavor singlet distribution, with q(x, µ 2 ) andq(x, µ 2 ) quark and antiquark distribution functions for flavor a; this contribution is absent in [9] and has been added by two of us in [35]. It was further found that a framework based on the Balitsky-Fadin-Kuraev-Lipatov [37][38][39] (BFKL) formalism and accounting both for quarks and gluons yields satisfactory agreement with data. However, there are remaining puzzles that are still left unanswered.
From a formal point of view, the above definition and its identification with the measured hadronic entropy has the obvious shortcoming that it relates an unphysical object, i.e. scheme dependent parton distribution functions (PDFs), to an observable, i.e. hadronic entropy. From the phenomenological side, the H1 collaboration determines hadronic entropy from the multiplicities of charged charged hadrons [40]. On the other hand, the theoretical framework of [9] is based on the treatment of purely gluonic emissions which give naturally rise to both charged and neutral hadrons. Comparing therefore the prediction of [9] with data for charged hadron multiplicities, one clearly expects a certain mismatch. A related issue is the large discrepancy between partonic and hadronic entropy encountered in [40], if the former is evaluated using leading order HERA PDFs.
The outline of this paper is as follows: in Sec. 2 we provide an overview of ambiguities and open questions in the current description and how they might reflect themselves in the comparison to data. In Sec. 3 we provide updated numerical results for the BFKL description of experimental results obtained by the H1 collaboration as well as a description based on leading order HERA PDF and PDFs obtained from an unintegrated gluon distribution subject to rcBK evolution. In Sec. 4 we provide a first analysis of the transition of this framework towards a region of phase space dominated by nonlinear QCD dynamics. In Sec. 5 we give our conclusions, while the appendix Sec. A summarizes details on the HSS and rcBK unintegrated gluon distribution.
2 Ambiguities in the current description

Overall normalization
The derivation of Eq. (1) is based on the expansion of the von Neumann entropy for large y = ln 1/x, or equivalently low x limit which is dominated by gluonic degrees of freedom. One finds S part. = − ∑ p n ln p n , p n (y) = e −∆y 1 − e −∆y n−1 , where p n denotes the probability to find n dipoles in the proton which satisfy p 1 (0) = 1 and p n>1 (0) = 0. They are obtained as a solution to the following evolution equation, with ∆ the BFKL intercept in the 2 dimensional model. Even within the 2-dimensional model, the above expression can be slightly generalized to with a certain constant C ≤ 1, which yields p n (0) ≤ 1 for all n. The mean number of dipoles is then obtained as Given that the gluon PDF is within the 2 dimensional model subject to the BFKL equation in zero dimensions, and taking into account that PDFs have a certain interpretation as number densities, i.e. their integral over momentum fraction x yields the expectation value of the parton number operator, it is somehow natural to interpret Eq. (5) as the gluon distribution. Even though xg(x) denotes usually the momentum fraction carried by gluons, while the number density is associated with the integral of g(x), the above identification is correct, since one really determines the mean value of the number of partons per ln(1/x), i.e. dn/dy , see also the more detailed discussion in Sec. 2.4.
There arises however an issue, whenever the identification of Eq. (5) as the PDF is lifted from the two dimensional model to four dimensions. While the PDF in four dimensions provides merely information about the one dimensional spatial structure of the proton, it carries an additional dependence on the factorization scale, which can be traced back to an integration over the two remaining transverse dimensions. It is this additional dependence on the factorization scale as well as the associated scheme dependence of PDFs which prohibit a direct identification of PDFs as number densities. Indeed, for renormalizable gauge theories the number density interpretation applies only to certain sum rules, i.e. the difference between the number of quarks and anti-quarks of a certain flavor, which is scheme independent, see e.g. [41] for a detailed discussion. For the gluon distribution such an interpretation may at best hold within e.g. light-cone gauge at leading order; beyond leading order, the PDF turns scheme dependent. In the light of such complication the best one can hope for is that Eq. (5) is proportional to the gluon distribution at leading order, with corresponding modifications, once higher order corrections are invoked. This observation leads us to the slightly modified relation, with some so far undetermined parameter B. Assuming 2 B =const. and lifting this expression to four dimensions, one would then realistically test in the comparison with data the evolution in y = ln 1/x and/or in ln Q 2 . Statements about the overall normalization should on the other hand be taken with care. Similarly the quark contribution is absent in the above expression due to the use of the purely gluonic effective 2 dimensional model. Since the low x seaquark distribution is driven by the gluon, the same statement applies to the relevance of this contribution, although it is naturally to be included as long as one interprets n as the mean number of partons. Sticking for the moment with the gluon distribution 3 , the above expression allows to re-express probabilities as which finally yields the following partonic entropy, Note that this expression is only meaningful for xg(x)/B > 1, which is the region where Eq. (8) yields 2 Within high energy factorization, B = const. is a meaningful assumption, since the x dependence is in general determined by the low x resummed gluon density; the same is true for collinear factorization at leading order. Collinear factorization beyond leading order (which is so far not worked out for this quantity) would most likely yield an x dependent parameter B; nevertheless, if the perturbative expansion is converging, this should imply a small correction. 3 To include quarks one would merely replace xg → xg + xΣ probabilities p n ≤ 1. Expanding Eq. (9) for xg(x) 1 and assuming B = O(1), one finds with e 2.71828 Euler's number. While for asymptotically large xg(x), the contributions due to B, e might be ignored, for realistic gluon distributions in the kinematic region covered by the HERA experiments, S part. ≤ 4 and while xg(x) 1 is satisfied, ln(xg) = O(1) and the above finite terms should be in principle kept.

Charged versus total hadron multiplicity
The above expression allows to determine partonic entropy as a function of the average number of partons in the system. It is then conjectured that the resulting partonic entropy agrees with the hadronic entropy. The latter can be obtained from the multiplicity of final state hadrons. Since the detection of neutral hadrons is experimentally challenging, this hadronic multiplicity distribution is usually determined for charged hadrons only. With pions the predominantly produced hadron species, and assuming that final state gluons turn with equal probabilities into positively, negatively, and neutral pion states, one can as a first estimate assume that the total number of produced hadrons is roughly 3/2 times the number of charged hadrons observed in experiment. In turn, the number of gluons and possibly seaquarks which yield charged hadrons is approximate the fraction 2/3 of the total parton number. This suggests to correct the partonic entropy by a corresponding factor, Clearly this factor is not exact, but merely an estimate of order of magnitude.

Gluon versus Quark contribution
In [9] entanglement entropy has been determined from a 2 dimensional model calculation, based on the dipole picture, which has been related to purely gluonic degrees of freedom. In a follow up study the same authors proposed in [36] to evaluate for the kinematic region covered by the HERA experiment the partonic entropy as the logarithm of the seaquark distribution. Finally, in [35], two of us proposed to evaluate partonic entropy as the logarithm of the sum of gluon and quark contribution. This treatment was motivated by the observation that the 2 dimensional gluonic model of [9] which identifies entropy as the logarithm of the average gluon number; for the complete theory, it seems therefore appropriate to generalize the average gluon number to the average number of quarks and gluons combined. While [35] managed to successfully describe HERA data, this description was plagued by several shortcomings. While the description based on the HSS gluon of [35] uses an inconsistent combination of overall normalization constants (which we correct in the subsequent numerical study, see also Sec. A for a detailed discussion), the description based on NNLO PDFs is subject to the above mentioned scheme dependence of parton distribution functions which allows for B = 1, as already point out in [35].
In the numerical study presented in Sec. 3, we show that a description based on a combination of quark and gluon contribution yields a good description of data within uncertainties if Eq. (11) is being employed, while we set B = e, i.e. we ignore all pre-asymptotic factors for the time being. While the observed agreement is pleasing, we believe that in the light of the above uncertainties in  12) for two different bin sizes. The result clearly depends on the size of the bin. The LO HERAPDF gluon distribution is evaluated at a factorization scale corresponding to the lower Q 2 value of the Q 2 bin as in [40] normalization, the main goal of the following study is the correct description of x and Q 2 dependence. The latter are directly related to evolution equations, to which the underlying distributions are subject to.

Binning and comparison to data
Assuming B = e for the moment, and assuming ln(xg) 1, partonic entropy is directly determined as the logarithm of the number of partons; indeed this has been the original the proposal of [9]. In practice this can however lead to confusion, since parton distribution functions are number densities, i.e. leaving aside the above mentioned issues of scheme dependence of such distributions, one expects that the number of gluons in a proton is obtained through with a similar expression in the case of quarks. In DIS reactions, the proton is on the other hand probed at a certain fixed value of Bjorken x. Experimentally this corresponds to counting final states for a certain bin size x ∈ [x min , x max ]. It is therefore tempting to define the number of gluons at a certain value of x through wherex ∈ [x min , x max ] and might be defined through x = x max x min dx xg(x, Q 2 ) x max or alternatively as the arithmetic mean of x min and x max . Partonic entropy is then obtained as ln n g (x).
To the best of our understanding, this is the method used by the H1 collaboration to determine the partonic entropy, with x max 5.71x min , see also the determination of this quantity using this method in Fig. 2. While the method yields the number of partons in the region of phase space [x min , x max ] if the PDF is interpreted as a number density, the result obviously depends strongly on the size of the interval [x min , x max ]. In particular the number of partons would approach zero, whenever the size of the interval turns infinitesimal small. For a meaningful comparison it is therefore necessary to study the number of partons, normalized to the bin size. Taking into account that usually one uses logarithmic bins in x, i.e. binning takes place in y = ln(1/x), one arrives at n g (x) = 1 y max − y min y max y min dy dn g dy = n g (y max ) − n g (y min ) y max − y min , y max,min = ln 1/x min,max , which in the limit of infinitesimal small bin sizes leads tō Despite of the usual association of xg(x, Q 2 ) with momentum sum rules, this quantity is therefore indeed the correct expression to compare to data binned in Q 2 and Y = ln(1/x). For the following numerical study, we further add the contribution due to quarks and average over the bin size in Q 2 , which finally yields the expression used for our comparison to data

Numerical results
Despite the above ambiguities in the overall normalization, we believe that is meaningful to compare at the current level of accuracy theory predictions to data. Ambiguities are mainly due overall normalization constants which for the partonic entropy turn into additive constants. In the following we use three theoretical models (to be described in more detail in the Appendix): • Leading order HERA PDF, subject to DGLAP evolution [42]. This particular set was chosen due to the observed mismatch between partonic and hadronic entropy by the H1 collaboration (which use LO HERAPDF). The goal of the comparison to LO HERAPDF is to demonstrate that such comparison should be done in the way outlined in the previous section, i.e. in bins of rapidity.
• Leading order PDFs calculated from the HSS unintegrated gluon, subject to BFKL evolution. The gluon density obtained in this scheme accounts for NLO corrections to the evolution kernel together with a collinear resummation of enhanced NLO contributions. The higher order corrections slow down the growth of the gluon density in the low x region, but do not lead to saturation. The fit is limited to the region Q 2 > 2 GeV 2 , see [31,32] for details.
• Leading order PDFs calculated from an unintegrated gluon, subject to rcBK evolution. The gluon density in this framework is subject to leading order BK evolution in x; albeit formally leading order, the evolution kernel includes NLO resummed running coupling corrections. The evolution takes into account effects due to high gluon densities and leads to a saturated gluon density (for developments that take into account exact kinematics in combination with saturation effects see [43][44][45]). A description of HERA data for the proton structure function F 2 is possible, since the nonlinear term in the evolution tames the rapid growth of the gluon distribution, induced by the linear term.
The first set has been used in [40] to compare to data. We show that once the corrective factor Eq. (11) is taken into account and the number of partons is evaluated as discussed in Sec. 2.4, this set of PDFs gives actually a good description of data, in contrast to the observation made in [40], leaving aside a small off-set in the normalization, which might be traced back to the effects discussed in Sec. 2.1.
The second set has been used previously in [35], while the last set is obtained from an unintegrated gluon distribution subject to Balitsky-Kovchegov (BK) evolution [46,47] which allows us to investigate possible contributions due to non-linear QCD evolution. In all three cases we use Eq. (1) to compare to data. The contribution due to seaquarks is in general small; the normalization of the description based on leading order HERA PDFs would slightly improve if one would only consider the contribution due to the gluon.
To obtain entropy from the low x QCD evolution equations, we need to determine PDFs i.e. integrated parton densities. The gluon density is obtained from where F (x, k k k 2 ) is the unintegrated gluon distribution, which is obtained from a solution to BFKL or BK evolution equations, including a fit to data. To obtain the quark PDF we apply the Catani-Hautmann procedure [48] xΣ(x, Q) where the splitting function reads [48] P qg z, k k k 2 and µ F denotes the factorization scale which we identify for the current study with the photon virtuality Q. k k k denotes the gluon transverse momentum and ∆ ∆ ∆ = q q q − zk k k with q q q the t-channel quark transverse momentum; In [35] the HSS integrated gluon Eq. (19) has been determined using a version of the unintegrated gluon with an overall normalization inconsistent with the normalization used for the determination of the quark PDF. In the following numerical study this has been corrected. In more general terms, the possibility of such inconsistencies can be traced back to a relatively large uncertainty in the overall normalization of the unintegrated gluon distributions, which have been obtained through fits which rely on the use of the leading order virtual photon impact factor. Since the latter is -as the seaquark distribution -proportional to α s , Results are compared to the final state hadron entropy derived from the charged multiplicity distributions measured by the H1 collaboration [40] for track pseudorapidities η * in the hadronic centreof-mass frame restricted to the range 0 < η * < 4.
it induces a relatively large normalization uncertainty on the extracted unintegrated gluon. To assess these uncertainties we multiply in the case of the HSS distributions the integrated gluon distribution by a factor α s (Q 2 )/α s (µ 2 ) and vary the renormalization scale µ in the range µ ∈ [Q/2, 2 · Q]. The unintegrated dipole gluon density subject to rcBK evolution is obtained through where N(r, x) is dipole amplitude obeying BK equation in coordinate space, see App. Sec. A for more details. In the above expression α s is kept constant, while it is running in the kernel of the evolution equation. In order to obtain uncertainties we use α s ∈ [0.2, 0.3], which are typical values for the hard scales investigated in this study. For the description based on leading order HERA PDFs we show the relatively small leading order PDF uncertainties, making use for the numerical evaluation of the package [49].
We find in Fig. 3 that the partonic entropy Eq. (18) obtained from the total number of partons gives a still a very good description of H1 data [40]. In particular within uncertainty bands, both HSS and rcBK give a good description of data. The HERA PDFs slightly overshoot the data but not as drastically as presented in [40]. This is due to the factor ln(2/3) which corrects for the fact that only charged hadrons have been measured as well as the effects discussed in Sec. 2.4. Since gluon and quark contribution add up, the description would slightly improve without the latter.

Towards the real photon limit
Having described the moderate and large Q 2 data, a natural question to ask is what happens if we go to lower values of Q 2 . With 1/Q 2 the effective area resolved in the interaction of virtual photon and proton, the limit Q 2 → 0 naturally leads to the case where the photon would observe the entire proton; entanglement entropy should be therefore absent in this limit. While the complete description of such a scenario is still unknown, we believe that investigating this limit within the current framework is already of interest. A more complete description should be probably formulated within the more general frameworks [50,51], where saturation effects [52] and therefore classicalization of wee gluons are taken into account [20]. For the moment we do not consider such modifications. Instead we will investigate the limit of small photon virtualities Q 2 . In Fig. 4 we show entanglement entropy as obtained from the solution to the rcBK evolution equation, as well as the HSS gluon and leading order HERA PDFs. While the rise with x flattens for three descriptions, this effect is clearly stronger for the rcBK description, which we attribute to effects related to gluon saturation. For leading order HERA PDFs and the HSS gluon, we limit the calculations to values Q 2 > 1 GeV 2 since for smaller values these descriptions would break down. This is however not a limitation for the rcBK or Golec-Biernat Wusthoff saturation model [53], which has been fitted for smaller values of Q 2 . From a formal point of view, one may justify this through the presence of a semi-hard saturation scale Q s (x) > Λ QCD . While results for Q 2 < 1 GeV 2 must be interpreted with care, they allow for a first qualitative investigation of this region of phase space. We will proceed with analytical studies of this region using the GBW model. Within this model, the dipole unintegrated gluon density reads 4 : 4 The BK for dipole gluon density has very similar x and k ⊥ dependence for |k ⊥ | < Q s where Q s = Q 2 0 x 0 x λ . The integrated gluon distribution is then obtained as 5 For large photon virtualities Q 2 Q s (x) 2 this yields i.e. the integrated gluon distribution is directly proportional to the saturation scale. In the limit of small photon virtualities one finds on the other hand: i.e. the integrated gluon distribution turns into a falling function of x. In particular one has the formal limit xg(x, Q 2 = 0) = 0. For entanglement entropy one finds, extrapolating the current framework to low Q 2 , Depending on the precise values of the transverse size S ⊥ , this expression will for some value of x turn eventually negative. Note however that the definition of probabilities Eq. (8) requires xg(x) ≥ B = e, for the current setup, which prevents us from reaching negative values. On the other hand for i.e. we recover the original expression for partonic entropy obtained in [9] plus a certain constant contribution.
A numerical study of the GBW model compared to rcBK is provided in Fig. 5. In the figure we plot only the contributions due to gluons. Entropy is being evaluated at low values of Q 2 < 1 GeV 2 since the both gluon densities have been fitted for this region of phase space. Moreover there exist in principle data for this region of phase and a comparison to those data would be of high interest, once they are in available in a suitable form. From a formal point of view, the existence of a semi-hard saturation scale allows for the evaluation of the gluon density at such low scales. We however stress here that the entropy formula in this region is an extrapolation. As can be seen from the plot to the right, the integrated gluon distribution is no longer necessarily large at lowest values of x. The expansion of Eq. (9) for large xg(x) is therefore not necessarily a good approximation. We therefore compare in Fig. 6 both the exact formula for entropy, derived within the 2 dimensional model, and its asymptotic expansion. The exact formula is not necessarily more accurate than the asymptotic expression, but the observed deviation indicates in which regions of phase space effects due to non-linear QCD dynamics might become relevant. In all of the the calculations we assume B = e which is the choice which minimizes the contribution due to constant terms. Note that the observation that for certain values of Q 2 the gluon PDF is a falling function of x can be directly linked to the behavior of the dipole amplitude featuring saturation stemming from a solution to BK or JIMWLK [54][55][56][57] evolution and which has been already observed for the solution of the BK equation in [58]. We interpret this as a mechanism of localization due to saturation of wee partons in longitudinal direction.

Conclusions
In this paper we continued the study of the proposal formulated in [9] to calculate for DIS reactions entanglement entropy from parton distribution functions. The purpose of this study has been twofold: on one hand we attempted to clarify some of the lose ends of this proposal such as the overall normalization, the relation between entanglement entropy and hadronic entropy of charged hadrons as well as a discussion on how to evaluate the number of partons for the determination of partonic entropy. On the other hand we provided a first exploration of the proposal towards low photon virtualities and very low x, where eventually non-linear QCD dynamics is expected to manifest itself.
To this end we provided a description of the data measured by the H1 collaboration through parton distribution functions subject to BFKL and BK evolutions equations as well as leading order HERA PDFs. As outlined in section 2, the description is at the moment mainly qualitative, due to various ambiguities in the theoretical description. In particular a precise phenomenology would require to work out an appropriate factorization of entanglement entropy into parton distribution functions and corresponding perturbative coefficients which would guarantee factorization scheme independence of the framework beyond leading order. While the frameworks used in this paper agree for moderate and large values of Q 2 , we find differences between BFKL evolution and leading order PDFs on the one hand and BK evolution on the other hand, if we turn to smaller values of Q 2 . In particular nonlinear dynamics, as encoded in the BK equation, suggests that entropy saturates or even decreases. On technical level we link this behavior to the feature of the BK gluon PDF which predicts a valence like gluon PDF i.e. a distribution falling as a function of x at some low enough value of x. This behavior of entanglement entropy is expected since at very low values of Q 2 the photon virtuality is so low that the system can not be bi-partitioned and entanglement entropy needs to vanish. To provide a definite answer to those questions, a generalization of the framework discussed here is necessary; in particular to explain the observed decrease of entropy with x. KK acknowledges the support of The Kosciuszko Foundation for the Academic year 22/23 for the project "Entropy of dense system of quarks and gluons". MH is grateful for hospitality at the Institute of Nuclear Physics and acknowledges support by Consejo Nacional de Ciencia y Tecnología grant number A1 S-43940 (CONACYT-SEP Ciencias Básicas).

A Some details on the rcBK and HSS
The procedure of getting the HSS unintegrated gluon density relies on solving NLO BFKL equation with BLM prescription for scale choice, collinear improvements and fitting initial condition to HERA data [31,32]. The rcBK is basically the LO BK equation formulated in the coordinate space and with running coupling constant included. The initial conditions are fitted [33] to HERA data [59].
• The unintegrated dipole gluon density from the BK is obtained via where N(r, x) is dipole amplitude obeying BK equation in the coordinate space. We use S ⊥ = 16.4 mb and α S is kept constant 6 . In order to assess normalization uncertainties due to the choice of this coupling constant, we vary it in the range [0.2-0.3]]. The BK equation reads [61] ∂N(r,Y ) ∂Y = d r 1 K run ( r, r 1 , r 2 )(N(r 1 ,Y ) + N(r 2 ,Y ) − N(r,Y ) − N(r 1 ,Y )N(r 2 ,Y )) where, Y = ln(x 0 /x) and r 2 = r − r 1 . The kernel K run is given by K run ( r, r 1 , r 2 ) = α s (r 2 ) 2π 2 N c r 2 r 2 1 r 2 2 + 1 r 2 1 α s (r 2 1 ) α s (r 2 2 ) with α s (r 2 ) = 12π the maximal allowed value for α s is α s,max = 0.7. The initial condition is given by McLerran-Venugopalan model [33,63] N MV (r, The equation was solved by newly developed code using Runge-Kutta method and using the parameters in the fit [33] to HERA data [59]. The numerical values of the parameter read Q s0 = 0.165GeV , γ = 1.135, C = 2.52, Λ QCD = 0.241, N f = 3, S ⊥ = 16.4 mb.
• The HSS gluon density follows from the BFKL evolution equation with the kernel obtained at NLL accurancy with collinear improvements [31,32,64]. The unintegrated gluon density reads where M is a characteristic hard scale of the process which we identify with Q.ĝ is an operator in γ space,ĝ x, whereᾱ s = α s N c /π with N c the number of colors and χ(γ) is the next-to-leading logarithmic (NLL) BFKL kernel which includes a resummation of both collinear enhanced terms as well as a resummation of large terms proportional to the first coefficient of the QCD beta function.
For the current study we set M = M = Q and n f = 4 with Λ QCD = 0.21 GeV. Q 0 = 0.28 GeV, and δ = 6.5. have been determined from a fit to the F 2 structure function in [31]. In this fit the overall running coupling constant has been evaluated at the renormalization scale µ 2 = QQ 0 , with Q the photon virtuality. For the construction of parton distribution function µ 2 = Q 2 is however a more natural choice. We therefore reevaluated the underlying fit and found that data on the proton structure F 2 [59] are equally well described, if we use µ 2 = Q 2 for the photon impact factor with a normalization C = 4.31. It is then this convention which we use in this study. In [35] the integrated quark distribution has been evaluated using C = 4.31, while the gluon has been evaluated with the normalization constant corresponding to the scale choice µ 2 = QQ 0 . Since distributions arise from the same unintegrated gluon distribution, such a treatment is not consistent. This has been corrected in the present study, including a variation of the overall running coupling to indicate the normalization uncertainty of our result.