Effective Debye Screening Mass in an Anisotropic Quark Gluon Plasma

Due to the rapid longitudinal expansion of the quark-gluon plasma created in heavy-ion collisions, large local-rest-frame momentum-space anisotropies are generated during the system's evolution. These momentum-space anisotropies complicate the modeling of heavy-quarkonium dynamics in the quark-gluon plasma due to the fact that the resulting inter-quark potentials are spatially anisotropic, requiring real-time solution of the 3D Schr\"odinger equation. Herein, we introduce a method for reducing anisotropic heavy-quark potentials to isotropic ones by introducing an effective screening mass that depends on the quantum numbers $l$ and $m$ of a given state. We demonstrate that, using the resulting effective Debye screening masses, one can solve a 1D Schr\"odinger equation and reproduce the full 3D results for the energies and binding energies of low-lying heavy-quarkonium bound states to relatively high accuracy. The resulting effective isotropic potential models could provide an efficient method for including momentum-anisotropy effects in open quantum system simulations of heavy-quarkonium dynamics in the quark-gluon plasma.


I. INTRODUCTION
The ongoing heavy-ion collision experiments at the Relativistic Heavy Ion Collider and the Large Hadron Collider aim to create and study a primordial state of matter called a quark-gluon plasma (QGP). One of the key observables that experimentalists are interested in is the production of heavy-quarkonium bound states in nucleus-nucleus (AA) collisions relative to their production in proton-proton (pp) collisions. This observable is of interest because it is expected that the formation of heavy-quarkonium bound states is suppressed in AA collisions relative to that in pp collisions. In early models the suppression of heavy-quarkonium bound states was solely due to the Debye screening of the inter-quark potential [1][2][3][4]. Since these early works, it has been shown that it is equally important to include effects of in-medium singlet-octet transitions and Landau damping of the exchanged gluons, which results in an imaginary contribution to the heavy-quark (HQ) potential [5][6][7][8][9][10][11][12][13][14]. These effects can be modelled systematically using numerical solutions of the Lindblad equation which governs the evolution of the in-medium heavy-quarkonium reduced density matrix [15][16][17][18].
One of the limitations of prior studies of the evolution of the heavy-quarkonium reduced density matrix is that they explicitly rely on an assumption of momentum-space isotropy. This assumption allows one to simplify the resulting dynamical evolution equations such that only solution of a one-dimensional Schrödinger equation (SE) subject to stochastic jumps is required [17]. This assumption can, however, not be made in a real-world quark-gluon plasma due to the rapid longitudinal expansion of the QGP and the existence of a finite relaxation time of the system. As a result, the QGP develops a high degree of momentum anisotropy in its local rest frame, see e.g. [19,20] and references therein. This complication means that, in practice, one must solve a 3D stochastic Schrödinger equation or corresponding Lindblad equation, which is numerically prohibitive. In this paper, we introduce a method to obtain effective 1D isotropic potentials from underlying anisotropic potentials. The method introduced can be used to simplify the inclusion of momentum-space anisotropy effects on heavy-quarkonium dynamics.
In order to develop and test this method, we will make use of a widely used anisotropic distribution function ansatz called the Romatschke-Strickland (RS) form [21,22] f aniso (k) ≡ f iso 1 λ where the isotropic distribution function f iso is an arbitrary isotropic distribution function that decreases sufficiently rapidly at large momentum. In the above equation, the unit vector n is used to denote the direction of anisotropy and λ is a temperature-like scale which, in the thermal equilibrium limit, should be understood as the temperature T of the system. In addition, we use an adjustable parameter ξ to quantify the degree of momentum-space anisotropy where k z ≡ k · n and k ⊥ ≡ k − n(k · n) correspond to the particle momenta along and perpendicular to the direction of anisotropy, respectively. By assuming n to be parallel to the beam-line direction, the case ξ > 0 is relevant to high-energy heavy-ion experiments. Among the many new phenomena which emerge as a consequence of QGP momentum-space anisotropy, herein we are interested in the in-medium properties of quarkonium states. In the non-relativistic limit, the binding energies and the decay widths of the bound states can be obtained by solving the 3D Schrödinger equation with a specified heavy-quark potential which describes the force between the quark and antiquark. Given the distribution function in Eq. (1), the complex-valued potential has been computed in the resummed hard-thermal-loop (HTL) perturbation theory [6,9,23]. 1 For arbitrary anisotropy ξ, the real part of the potential has to be evaluated numerically and analytical results can only be obtained in the limits, e.g., |ξ| ≪ 1, ξ ≫ 1, and ξ → −1. On the other hand, the determination of the imaginary part for general ξ has encountered difficulties related to the presence of a pinch singularity in the static gluon propagator [13]. As a result, this is still an open question and we will focus on the real part of the potential in the current work.
It is worth noting that the perturbative HQ potential in an anisotropic plasma can be well fitted with a standard Debye-screened function ∼ e −rmD(λ,ξ,θ) /r, therefore, is formally analogous to its counterpart in equilibrium [24]. The determination ofm D (λ, ξ, θ) requires a specific angle θ denoting the quark pair alignment r with respect to the direction of anisotropy n. However, such a θ-dependence results in an unclear physical interpretation of the screening mass and only makes sense in a classical picture where the quark pair has a fixed alignment. In addition, solving the three-dimensional SE makes the numerical determination of the binding energies of quarkonium states much more complicated compared to the case where a spherically symmetric HQ potential can be used. At present, this is the main obstacle to developing phenomenological applications which include momentum-anisotropy effects.
In this work, we propose an angle-averaged screening mass M lm (λ, ξ) in an anisotropic medium through the following definition where Y lm are spherical harmonics with the azimuthal quantum number l and magnetic quantum number m. Notice that although no angular dependence appears in M lm (λ, ξ), this quantity is not universal for all the states labelled by the set of quantum numbers l and m. As will be discussed in this work, information on the physical properties of quarkonium states in an anisotropic QCD plasma can be obtained at a quantitative level by analyzing the corresponding problem in an "isotropic" medium characterized by the angle-averaged screening mass M lm (λ, ξ). In this sense, M lm (λ, ξ) is also called the effective screening mass. The rest of the paper is organized as the following. In Sec. II, we introduce the model construction of the HQ potential in an anisotropic QCD plasma and qualitatively explain the rationality of our definition of the effective screening mass as given by Eq. (3). In Sec. III, an anisotropic HQ potential as given in Eq. (4) is Taylor expanded around an isotropic configurationm D (λ, ξ, θ) = M lm (λ, ξ) where the values of l and m are specified based on the concept of the "most similar state". In Sec. IV, we demonstrate that highly accurate results of the eigen/binding energies of quarkonium states can be obtained by solving the Schrödinger equation with only the leading order contribution in the above Taylor expansion of the HQ potential. In addition, an estimation on the discrepancy resulting from replacingm D (λ, ξ, θ) with M lm (λ, ξ) in Eq. (4) is also obtained. Some applications of the effective screening mass derived herein are considered in Sec. V. We briefly summarize our findings and give an outlook for the future in Sec. VI. Finally, for several low-lying heavy-quarkonium bound states, the exact 3D results of the eigen/binding energies based on the potential model in Eq. (4), together with the corresponding discrepancies from using the one-dimensional potential model based on the effective screening masses are listed in Appendix A, which provides a direct numerical check of our method.

II. THE HEAVY-QUARK POTENTIAL IN AN ANISOTROPIC MEDIUM WITH EFFECTIVE SCREENING MASS
Previous works have demonstrated that the non-perturbative contributions to the HQ potential are not negligible for practical applications, therefore, the potential derived from perturbation theory can not fully capture the inmedium properties of quarkonium states. In general, describing the long-distance interaction between the quark pair relies on the phenomenological potential models which are constructed based on the lattice simulations. As a popular research topic, continuous attention has been paid on modeling the real part of the potential. In addition, developing a complex potential model has also been considered in recent years. However, most of the studies on quarkonium physics were limited to a momentum-space isotropic quark-gluon plasma. Due to the absence of the lattice simulations for a non-ideal anisotropic system, the corresponding potential models have to be introduced in an indirect way.
A first attempt to model the heavy-quark potential in an anisotropic plasma was carried out in Ref. [10]. As a "minimal" extension of the isotropic Karsch-Mehr-Satz (KMS) potential [1,2] based on the internal energy, the anisotropic version was obtained by replacing the ideal Debye mass m D (λ) in thermal equilibrium with an anisotropic screening scalem D (λ, ξ, θ). Explicitly, one assumes the following form for the (real part of the) potential model where α is an effective Coulomb coupling at (moderately) short distance and σ is the string tension. In the rest of the paper, we choose α = C F α s = 0.385 and σ = 0.223 GeV 2 for numerical evaluations. These two parameters are assumed to be unchanged in a hot medium, once determined at zero temperature. Following Ref. [10], the isotropic Debye mass is given by m D = Agλ 1 + N f /6 with N f = 2 and g = 1.72. The factor A = 1.4 which accounts for non-perturbative effects has been determined in lattice calculations. In addition, we assume the critical temperature T c = 192 MeV. For our purposes the most important feature of the KMS potential model is that medium effects are entirely encoded in the Debye mass m D (λ) and the very same screening scale that appears in the perturbative Coulombic term also shows up in the non-perturbative string contributions. Extending this assumption to anisotropic plasmas, the feasibility of the "minimal" extension depends on the extraction of an angle-dependent screening massm D (λ, ξ, θ) from the perturbative potential which has been computed from the first principles. Despite this choice of potential, we emphasize that the generalization from an isotropic HQ potential to the anisotropic one is not restricted to the use of the KMS potential model. In fact, any other isotropic model, as long it shares this feature with the KMS model, can be generalized in a similar way. Some recent examples of such HQ potential models can be found in Refs. [14,[25][26][27][28]. For definiteness, in the following discussions, we will consider the potential model as given in Eq. (4).
From a phenomenological point of view, we expect that non-zero anisotropy corrections only amount to a modification on the isotropic Debye mass m D (λ). In fact, this turns to be reasonable based on an analysis of the perturbative contributions in the HQ potential. In principle, the anisotropic screening massm D (λ, ξ, θ) can be defined as the inverse of the distance scale over which |rV (r)| drops by a factor of e. For small anisotropy ξ, the screening mass is given by [10]m As expected, the non-ideal Debye massm D (λ, ξ, θ) depends not only on the anisotropy parameter ξ but also on the angle θ between r and n. With this in hand, we can approximate the anisotropic potential at short distances with a Debye-screened Coulomb potential where the θ-dependent screening mass is given by Eq. (5). More interestingly, it is found that for arbitrary anisotropies, the perturbative potential can also be perfectly parameterized by using the same Debye-screened Coulomb potential although in this case, the anisotropic Debye massm D (λ, ξ, θ) takes a much more complicated form as [24]m where The above expression ensures the correct asymptotic limits for small and large ξ and efficiently reproduces the shortrange anisotropic potential when compared to the direct numerical evaluation of the full result [24]. The existence of such a parameterization further guarantees the justification of the model construction in an anisotropic medium as proposed above. Namely, the behaviors of the potential at short distances can be well described by the Debye-screened Coulomb potential from which an anisotropic screening mass can be extracted.
Although it is not surprising that a θ-dependence emerges in the anisotropic screening mass, one has to face the technical difficulty of solving a two-or three-dimensional SE. On the other hand, if we replacem D (λ, ξ, θ) with the effective screening mass as introduced in Eq. (3), the restoration of the spherical symmetry in the HQ potential greatly simplifies the numerical evaluations as one only needs to solve a one-dimensional problem. We believe this is very important for many studies concerning quarkonium physics in an anisotropic situation. Therefore, the key issue is to demonstrate the rationality of the definition of M lm (λ, ξ) introduced in Eq. (3). Qualitatively, it can be understood by analyzing both the anisotropy and quantum number dependence of the effective screening mass.
Using Eq. (5), for the anisotropic screening mass in small ξ limit, the integrations in Eq. (3) can be performed analytically, leading to with As expected, the effective screening mass decreases with increasing anisotropy ξ and a reduced screening effect exists in an anisotropic plasma since the values of M lm (λ, ξ) are always smaller than the ideal Debye mass m D (λ). Furthermore, the observed polarization of quarkonium states with non-zero angular momentum in an anisotropic plasma can be qualitatively explained by the corresponding m-dependence as given in Eq. (8). We find that instead of m D (λ)(1 − ξ/5) for P 0 states, M lm (λ, ξ) takes a different value m D (λ)(1 − 3ξ/20) for P ±1 states. As a result, one can expect that the first excited states of bottomonia 2 with L z = 0 have a smaller screening mass therefore is bounded more tightly as compared to those with L z = ±1. For arbitrary ξ, the above discussions also hold although in this case, the values of M lm (λ, ξ) need to be determined numerically. In addition, since the function k(l, m) is in the range of 1 ≤ k(l, m) ≤ 8/5, Eq. (8) suggests that at a given temperature-like scale λ and anisotropy ξ, the effective screening mass has a minimum value m D (λ)(1 − ξ/5) with l = 1 and m = 0 while its maximum, equaling to m D (λ)(1 − ξ/8), is reached when l → ∞ and m = ±l. Numerical evaluations on the values of M lm (λ, ξ) for arbitrary ξ turn to support the same conclusion and we give the corresponding results with anisotropies ξ = 5 and ξ = 10 for different l and m in Table I and Table II, respectively.

III. PERTURBATIVE EVALUATIONS ON THE ENERGIES OF QUARKONIUM STATES
The above discussions suggest that it is reasonable to introduce an angle-averaged screening mass as defined in Eq. (3) to replace the angle-dependent mass in Eq. (4). However, the discrepancy resulting from such a replacement has to be estimated at a quantitative level to ensure an accurate determination on the physical properties of quarkonium states in an anisotropic plasma. We start with the stationary Schördinger equation where m Q is the reduced mass for the quarkonium bound state with eigen-energy E andL 2 is the square of the angular-momentum operator. For an isotropic potential V (r) without the dependence on the azimuthal angle θ and polar angle φ, the wave-function ψ(r) can be separated into a radial and an angular part as and the associated eigen-energy denoted by E nl is determined by the following equation In general, Eq. (12) can not be solved analytically but the numerical evaluations can be simply carried out. According to the well-known nodal theorem, the level of excitations can be identified by the number n of the nodes of the reduced radial wave-function R nl (r). For a given l, the ground-state wave-function corresponds to n = 0, while the wavefunction of some radial excitations has a non-zero and finite number of nodes. Therefore, the number n = 0, 1, 2, 3 · · · . In addition, the spherical harmonics Y lm (θ, φ) in Eq. (11) is standard with m = −l, −l + 1, · · · , 0, · · · , l − 1, l and l = 0, 1, 2, 3, · · · . Because of the breaking of the spherical symmetry, for an anisotropic potential V (r), the above separation of the wave-function is no longer true, namely, ψ iso nlm (r) is not an eigen-state for the system. On the other hand, due to the completeness of the eigen-functions in Eq. (11), we can expand ψ aniso (r) as the following where k is used to label the level of the excitations. Accordingly, the eigen-energy associated with ψ With the effective screening mass proposed in Eq. (3), the anisotropic potential as given in Eq. (4) can be expanded around an isotropic configurationm D (λ, ξ, θ) = M l ′ m ′ (λ, ξ) with the leading order contribution is denoted by . The expansion can be formally expressed as the following Taylor series where the radial and angular dependences in the anisotropic terms can be separated as with If the Taylor series converges quickly, higher order terms V (s) l ′ m ′ (r) for s ≥ 1 act as a small perturbation to the isotropic potential V (0) l ′ m ′ (r). As a result, the eigen-energy E nl determined by Eq. (12) with V (r) = V (0) l ′ m ′ (r) could be an ideal approximation to the eigen-energy E [k] of the bound state in an anisotropic medium. However, there exists an ambiguity due to fact that both the reduced radial wave-function R nl (r) and the eigen-energy E nl are not unique because the above isotropic potential constructed with the effective screening mass has a (l ′ , m ′ )-dependence. For example, with different values of l ′ and m ′ , a set of ground states ψ iso 000 (r) can be obtained based on Eqs. (11) and (12). A question naturally arising is which values of l ′ and m ′ lead to an eigen-energy E 00 that is closest to E [0] . This question can be answered by looking at the energy correction ∆E which in the first approximation in the perturbation theory reads Naively, one can expect that the main contribution to the energy correction comes from the first term ∆E 1 . Therefore, when choosing (l ′ , m ′ ) = (0, 0), this term vanishes by definition and the perturbative correction to the eigen-energy E 00 becomes negligible. The above conjecture indicates that the isotropic potential V (0) 00 (r) constructed with the effective screening mass M 00 is the right one that should be used in Eq. (12). The resulting ground state described by the wave-function ψ iso 000 (r) is the "most similar state" in the sense that the associated eigen-energy E 00 is the best approximation to the lowest energy E [0] in an anisotropic medium. Furthermore, using the "most similar state" is actually consistent with the explanation for the polarization of states with non-zero angular momentum appearing in an anisotropic plasma which has already been discussed in Sec. II. This conclusion can be generalized to excited states. Solving the stationary SE with an isotropic potential V (0) l ′ m ′ (r), among the eigen-states as formally given by Eq. (11), we are interested in a set of the states with the quantum numbers (l, m) = (l ′ , m ′ ). The eigen-energy E [k] can be well approximated by E nl ′ of the so-called "most similar state" ψ iso nl ′ m ′ (r). In Table III, we list the eigen-energies as well as the perturbative corrections evaluated with different M l ′ m ′ for the bottomonium states (including Υ, χ b0 and χ b±1 ) at ξ = 1 and λ = 1.1T c . Comparing with the exact values, .318 MeV and E [2] = 507.878 MeV which are obtained by solving the threedimensional SE with the anisotropic HQ potential in Eq. (4), we find that the energies E 00 or E 01 associated with the "most similar state" are closest to exact results as expected. The above discussion shows that one-dimensional SE should be sufficient to provide the desired information on the bound states in an anisotropic QCD medium at a quantitative level, namely, the solution with an isotropic potential V (0) l ′ m ′ (r) is indeed a very good approximation to that of an anisotropic medium characterized by a screening massm D (λ, ξ, θ). In addition, finding excited states in a numerical approach is also a challenging task especially when the anisotropic situation is taken into account. On the other hand, when solving the Schrödinger equation with the (l ′ , m ′ )-dependent potential V (0) l ′ m ′ (r) , the obtained "ground state" ψ 0l ′ m ′ (r) for (l ′ , m ′ ) = (0, 0) actually corresponds an excited state in the anisotropic medium. For example, the first two excited states of bottomonia denoted as χ b0 and χ b±1 are identified with the "ground states" determined by V (0) 10 (r) and V (0) 1±1 (r), respectively. Therefore, finding the excited states can be also simplified in our approach based on the effective screening mass.
115  Furthermore, the Taylor series given in Eq. (14) is very different from the expansion of the anisotropic potential V (r) around ξ = 0. Although the leading contributions are both independent of the angles, the energy splitting of states with angular quantum number l = 0 already shows up at leading order due to the m-dependence of the effective screening mass when Eq. (14) is considered. On the other hand, in order to see such a splitting of energy in the latter approach, one has to take into account the anisotropic higher order corrections. Of course, the main disadvantage of the expansion around ξ = 0 is the limitation for application at large anisotropies.

IV. DISCUSSION ON THE SUPPRESSION OF THE TAYLOR SERIES
The "most similar state" ψ iso nl ′ m ′ (r) is an eigen-state determined by Eqs. (11) and (12) where an isotropic heavyquark potential V (r, M l ′ m ′ ) should be taken into account. As shown in Table III, the perturbative corrections to the eigen-energy are so small that ψ iso nl ′ m ′ (r) could be an ideal approximation to the corresponding eigen-state ψ [k] aniso (r) in an anisotropic plasma. Therefore, the one-dimensional SE with the effective screening mass can be used to study quarkonia in an anisotropic medium. For any given quantum numbers n, l and m, the eigen-energy corrections can be written as Notice that different from Eq. (17), the sum in the above equation starts from s = 2 because only the "most similar state" will be considered in the following. In order to understand the reason why the contribution from Eq. (18) is negligible, we formaly divided ∆E into the following two parts where and A. The angular part of the perturbative corrections to the eigen-energies We first look at the angular part. Since the spherical harmonics are universal, Eq. (20) can be determined once the exact form of the anisotropic Debye mass is known. This means ∆E ag s is independent of the specific potential model used in the Schrödinger equation. Starting with the case where ξ is small, we have where x ≡ cos θ. One can expect that with increasing s, the (absolute) values of ∆E ag s (ξ ≪ 1) get decreased simply due to the factor (ξ/8) s . In fact, we can further show that there is an extra suppression from the (l, m)−dependent term in the above equation which is denoted as Here, C n s is the binomial coefficient and X (n) is defined as X (n) = Y lm |x n |Y lm . Using the iterative formula for the spherical harmonics it is straightforward to obtain an analytical expression for X (n) , although this is rather tedious for large n. Explicitly, for s = 2, we arrive at with X (2) = a 2 lm + a 2 l−1,m and X (4) = a 2 lm (a 2 l+1,m + a 2 lm + a 2 l−1,m ) + a 2 l−1,m (a 2 lm + a 2 l−1,m + a 2 l−2,m ) .
In fact, one can prove that Eq. (25) approaches zero when l → ∞ and m = ±l and its maximum value 68/441 can be obtained when l = 2 and m = 0. Therefore, I lm ≤ 68/441 ≈ 0.154. For large s, the (l, m)−dependence in I (s) lm turns to be very complicated. Instead, we will consider the following ratio for n = 1, 2, 3, · · · |I (2n+1) lm Here, δ lm is defined as the maximum value of |x 2 + 1 − k(l, m)| for a given (l, m) when x 2 changes from 0 to 1. Since 1 ≤ k(l, m) ≤ 8/5, δ lm can not be larger than 1. In fact, the largest δ lm is given by δ ∞,±l = 1. In addition, lm < δ 2 lm can be also shown similarly. Based on the above discussion, we conclude that the magnitude of I (s) lm is small and gets suppressed when s is large. In particular, the maximum value of I (2) lm is given by 0.154. For s > 2, although the exact values of the maximum are not computed here, it is possible to estimate an upper limit which is given by |I We should mention that for even s, I lm is monotonically decreasing with increasing s, however, the same conclusion can not be drawn for s = 2, 3, 4, · · · . Instead, only I (2n) lm > |I (2n+1) lm | for n = 1, 2, 3, · · · can be justified. Furthermore, the iterative formula in Eq. (24) turns to be very useful to analytically study the angular part of the perturbative energy correction even for the case of arbitrary anisotropies as long as the anisotropic Debye mass can be parameterized as a polynomial function of x, namelym D (λ, ξ, θ) = n i=0 c i (λ, ξ)x i . However,m D (λ, ξ, θ) as given in Eq. (6) doesn't satisfy this condition.
For arbitrary ξ, we start with the simplest case where s = 2 In the above equation f (ξ) = f 2 (ξ)/f 1 (ξ) and the explicit expression of f (ξ) is rather complicated which we don't list here. However, the key point is that f (ξ) is positive for ξ > 0 and it is not a monotonic function of ξ. The minimum can be numerically determined as f min (ξ) ≈ f (2.958) ≈ 0.694. In Fig. 1, we plot f (ξ) as a function of ξ. It is not easy to get a general expression of ∆E ag 2 (ξ) for any given l and m. On the other hand, in order to see the non-trivial ξ-dependence of ∆E ag 2 (ξ), we fix l and m. The simplest case is to consider (l, m) = (0, 0) where the numerator in Eq. (29) can be carried out analytically and the result takes the following form where the integral in the denominator corresponds to a hypergeometric function. Numerically, it can be shown that the above result decreases with increasing values of f (ξ), therefore, the maximum of ∆E ag 2 (ξ) with (l, m) = (0, 0) appears when ξ ≈ 2.958, see Fig. 2 (the left panel). Changing the values of l and m, we have checked that the maximum of ∆E ag 2 (ξ) is always at ξ ≈ 2.958 and the maximum values of ∆E ag 2 (ξ) at a given (l, m) are list in Table IV. Our numerical results also suggest that as a function of l and m, ∆E ag 2 (2.958) becomes largest when (l, m) = (2, 0) which coincides with the finding in the small ξ limit.   For s > 2, we are not going to evaluate ∆E ag s (ξ) with various quantum numbers (l, m), instead, we try to estimate the corresponding upper limit which should be sufficient to demonstrate the suppression of ∆E ag s (ξ) when s is large. It is clear that the function is monotonically decreasing with x 2 when x changes from 0 to 1, therefore, the maximum of F 2 lm (ξ, x) locates at either x = 0 or x = 1 which can be denoted as max(F 2 lm (ξ, x = 0), F 2 lm (ξ, x = 1)). Therefore, the ratio of ∆E ag 2n+2 (ξ)/∆E ag 2n (ξ) for n = 1, 2, 3, · · · , satisfies Then, the question amounts to estimating an upper limit for F 2 lm (ξ, x = 0) and F 2 lm (ξ, x = 1) for arbitrary anisotropies and the quantum numbers. Notice that f (ξ) > 0 and f min (ξ) ≈ 0.694, we have The estimated upper limit is valid independent on the values of (l, m). The analysis on F 2 lm (ξ, x = 1) can be carried out in a similar approach and we can show On the other hand, once l and m are specified, F 2 lm (ξ, x = 0) and F 2 lm (ξ, x = 1) can be evaluated numerically and one can further reduce this limit. As shown in Fig. 2 (the right panel), with (l, m) = (0, 0), the maxima of F 2 lm (ξ, x = 0) and F 2 lm (ξ, x = 1) which appear at ξ = 2.958 are found to be 0.0080 and 0.0165, respectively. More results are given in Table V. At this point, the estimation that max(F 2 lm (ξ, x = 0), F 2 lm (ξ, x = 1)) < 0.0625 is justified. In addition, the square root of this upper limit corresponds to that of the ratio |∆E ag 2n+1 (ξ)|/∆E ag 2n (ξ), namely, the following relations are always true regardless of the specific values of ξ as well as the quantum numbers l and m, ∆E ag 2n+2 (ξ) < 0.0625∆E ag 2n (ξ) , and |∆E ag 2n+1 (ξ)| < 0.250∆E ag 2n (ξ) .
For practical applications, we are probably more interested in quarkonium states with not very large azimuthal quantum number. Based on the results given in Tables IV and V, we can show the following The discussions above of the angular part of the perturbative corrections depend only on the parameterization of the anisotropic Debye massm D . On the other hand, one has to specify an explicit form of the heavy-quark potential in order to study the corresponding radial part ∆E ra s as defined in Eq. (21). In addition, unlike the spherical harmonics Y lm (θ, φ), the radial wave-function R nl (r) can, in general, only be obtained numerically by solving Eq. (12) .
To proceed further, we consider Eq. (4) as the heavy-quark potential model. The explicit form of G (s) lm defined in Eq. (16) can be written as To obtain the above equation, we have used the following derivatives In Fig. 3, we show G (s) lm as a function ofr ≡ rM lm for s = 2, 3, 4, 5. The plot on the left corresponds to the effective screening mass M lm = 1500 MeV and we focus on a region of the dimensionless variable up tor = 3.75. Therefore, the size of the quarkonium states are assumed to be smaller than 0.5 fm, which roughly speaking is the upper limit of the root-mean-square radii of quarkonia above the critical temperature. The right plot shows the results at M lm = 500 MeV. Correspondingly, we consider a relatively narrow region ofr. For the "most similar state", due to the vanishing angular part in the perturbative correction, ∆E 1 = 0 which is independent of the values of G (1) lm . On the other hand, the plots in Fig. 3 suggest that the magnitude of G (2) lm doesn't exceed ∼ 70 MeV within the given region ofr that is relevant to quarkonium studies 3 . Therefore, ∆E 2 is estimated to be less than 0.6 MeV given that ∆E ag 2 ≤ 0.00823, see Table IV. In the same region ofr, the magnitude of G (s) lm gets smaller as s is increasing, naively, we don't expect any enhancement from the radial part to ∆E s when s is large.
As a very rough estimation, the above discussion presumes that R nl (r) is simply a Dirac delta and becomes nonvanishing only if r equals the root-mean-square radius of the bound state. To do it in a relatively rigorous way, we consider a normalized function This is actually the radial part of the wave-function for ground state when the potential is Coulombic. However, the Bohr radius a 0 should be understood as the most probable radius of a quarkonium state. For the Coulomb term, one obtains ∆E ra,C withâ 0 ≡ M lm a 0 . Here, we have used the following integral where a > 0 and n is a positive integer. Similarly, the contribution from the string term is given by Notice that after integrating over r with the help of Eq. (41), the summation over n appearing in Eq. (37) can be carried out by using the identity s n=0 (n + 1)(n + 2) Both ∆E ra,C s and ∆E ra,S s have non-trivial dependences on the most probable radius a 0 as well as the effective screening mass M lm . In order to estimate an upper bound for the the radial part of the perturbative corrections |∆E ra s | ≡ |∆E ra,C s + ∆E ra,S s |, the value ranges of a 0 and M lm have to be specified. When the screening effect becomes strong enough, even the smallest quarkonium state Υ can not survive. Therefore, it is reasonable to assume M lm ≤ 1500 MeV which, based on the potential model adopted here, corresponds to a temperature less than 3T c in the limit ξ = 0. On the other hand, the most probable radius a 0 equals r 2 / √ 3, as a result, a 0 = 0.3 fm is equivalent to a root-mean-square radius r 2 ≈ 0.5 fm which has been considered as the largest size of quarkonia that could be bound in the deconfining phase.
As demonstrated in the left plot of Fig. 4, the magnitude of ∆E ra 2 is always smaller than ∼ 40 MeV when we vary the radius at a given M lm . Numerically, we find that the maximum of |∆E ra 2 | appears as both variables equal the largest values that they may take, namely, a 0 = 0.3 fm and M lm = 1500 MeV. As compared with our previous analysis, the upper bound of ∆E 2 now is reduced to ∼ 0.3 MeV. Although small enough in magnitude, we would like to take it as an overestimation. In fact, for any specific quarkonium state, in order to get ∆E ra 2 at any fixed M lm , the corresponding root-mean-square radius or a 0 has to be determined by (numerically) solving the Schrödinger equation with an isotropic HQ potential V (r, M lm ). According to the obtained values of the most probable radius a 0 , one can locate a point on each curve in Fig. 4 which we refer to as the "physical point". Taking the Υ state as an example, the "physical point" which corresponds to the physical a 0 of the quarkonium state in consideration, is denoted by a filled circle in Fig. 4.
In addition, under the constraint conditions that a 0 ≤ 0.3 fm and M lm ≤ 1500 MeV, the maximum of |∆E ra s | can be determined for any given s. In general, the maximum appears when both a 0 and M lm take their largest possible values. However, there are also exceptions such as |∆E ra 5 | max = |∆E ra 5 (a 0 ≈ 0.124 fm, M lm = 1500 MeV)|, see the right plot of Fig. 4. The corresponding results as presented in Fig. 5 suggest that as compared to |∆E ra 2 | max , higher order terms ∆E ra s>2 rapidly decrease in magnitude and become negligibly small even at moderate s. In fact, both Eq. (40) and Eq. (42) vanish when s is infinitely large. Notice that the sequential suppression of |∆E ra s | can not be guaranteed because the contribution from ∆E ra 2 can be zero as long as certain values are assigned to a 0 and M lm , see Fig. 4. On the other hand, due to the strong suppression on |∆E ra s | max as well as the angular part |∆E ag s |, the perturbative corrections ∆E as defined in Eq. (18) can be neglected within an error about ∼ 0.3 MeV.
According to Table III, although the differences among the eigen-energies of bottomonium obtained with various effective screening mass are not significant, the necessity of using the "most similar state" could be understood more clearly when we look at the energy splitting of quarkonium states with non-zero angular momentum which is a unique feature arising in an anisotropic medium. Without choosing the "most similar state", the contribution from perturbative correction could be on the order of ∼ 1 MeV which is larger than the upper bound ∼ 0.3 MeV estimated for the "most similar state" and actually comparable to the splitting between the states with different polarizations. For example, solving the exact three-dimensional SE, we find that the eigen-energy of χ b0 differs from that of χ b±1 by an amount ∼ 1 MeV. Therefore, the "most similar state", which is expected to provide the most accurate results, is essential to reducing the three-dimensional problem effectively to an isotropic one-dimensional problem. Excellent agreement can be seen in Tables VII, VIII, and IX when comparing the corresponding eigen-energies with the exact solutions.

C. The heavy-quark potential at large distances
For quarkonium studies, the binding energies of the bound states are of particular interest. Therefore, the behavior of the heavy-quark potential at large distances is essentially important since the binding energy, by definition, equals the eigen-energy minus V ∞ (λ, ξ). In an anisotropic plasma, V ∞ (λ, ξ) is given by the following expectation values where in principle the wave function ψ aniso (r) has to be determined by solving the three-dimensional SE. On the other hand, expanding V (r → ∞,m D (λ, ξ, θ)) around the effective screening mass, we find that the leading term V (0) ∞ (λ, ξ) which is independent on the azimuthal angle θ turns to be an ideal approximation of V ∞ (λ, ξ) when the "most similar state" is considered.
Using Eq. (4) as the heavy-quark potential model, V ∞ (λ, ξ) can be expressed as the following Taylor series In the first line of the above equation, we have used the "most similar state" denoted as ψ nlm (r) to approximate the corresponding wave-function ψ [k] aniso (r) in an anisotropic plasma. As a result, the s = 1 term in the above sum is absent. The higher order corrections with s ≥ 2 are suppressed as compared to the leading contribution V (0) ∞ (λ, ξ) = 2σ/M lm (λ, ξ) by a small factor of ∆E ag s . In fact, we can further show the following result according to Eq. (36) where . Because the heavy-quark potential at finite temperature is deeper than the vacuum potential, one can expect that V ∞ (λ, ξ) and V (0) ∞ (λ, ξ) are smaller than V ∞ (λ = 0). Roughly speaking, V ∞ (λ = 0) corresponds to the value of the vacuum potential at the string breaking distance which, in general, is approximately 1 GeV. Therefore, one can estimate that the higher order correction ∆V ∞ (λ, ξ) < 10 MeV.
894 .   In addition, the necessity of using the "most similar state" could be further justified when we expand V (r → ∞,m D (λ, θ, ξ)) around the effective screening mass M l ′ m ′ (λ, ξ) with different quantum numbers l ′ and m ′ . The leading term V (0) ∞ (λ, ξ) = 2σ/M l ′ m ′ (λ, ξ) while the corrections can be denoted as Taking (l, m) = (0, 0) as an example, the Taylor series of V ∞ (λ, ξ) has been computed and the results are shown in Table VI. The "most similar state" corresponds to using the effective screening mass M 00 (λ, ξ). In this case, V ∞ (λ, ξ) turns out to be a good approximation to V ∞ (λ, ξ) which is 896.20 MeV from the exact 3D evaluation. On the other hand, with other effective screening mass determined with (l ′ m ′ ) = (0, 0), the non-vanishing V (1) ∞ (λ, ξ) may give rise to a considerable correction on the binding energy, especially for excited states. More importantly, based on the estimation in Ref. [10], the splitting of the binding energy of χ b with different angular momentum is on the order of 50 MeV which is obviously comparable with the contribution from V (1) ∞ (λ, ξ). As a result, when keeping only the leading term in Eq. (45), one has to make use of the "most similar state" in order to get a quantitatively reliable result.
V. SOME APPLICATIONS Perturbatively, the Debye mass m D at leading order is proportional to the temperature λ or T in an isotropic (equilibrium) medium. On the other hand, studies on quarkonium states become most relevant in a temperature region not from far above T c where non-perturbative physics plays an important role. As argued in Ref. [29], lattice simulations suggest that the non-perturbative Debye mass can be approximated as a constant factor A times the perturbative m D , therefore, as used in our numerical evaluations, m D is also proportional to λ at relatively lower temperatures, roughly speaking, from T c to ∼ 3T c , which is known as the semi-QGP [30][31][32][33]. With the effective screening mass M lm (λ, ξ) given in Eq. (3), it is possible to define an effective temperatureλ in an anisotropic medium by which the effective screening mass can be formally expressed as M lm (λ, ξ) = Agλ 1 + N f /6. The effective temperatureλ is determined by requiring that the screening mass M lm (λ, ξ) or the binding energy of a bound state in an anisotropic medium is equal to that in an isotropic medium determined at temperatureλ. Apparently,λ depends not only on the anisotropy ξ but also on the quantum numbers l and m. According to Eqs. (6) and (3), we can show thatλ The ratioλ/λ as a function of ξ at some fixed (l, m) is shown in Fig. 6. As we can see, a linear decrease with the increasing anisotropies appears in the small ξ region. This is actually in accordance with Eq. (8) by which the above ratio reduces to 1 − k(l, m)ξ/8 for ξ ≪ 1. The perturbative Debye mass at non-zero quark chemical potential µ is given by withμ ≡ µ/(2πλ), which suggests an enhanced screening in the high temperature limit where the potential takes a Debye-screened form. Therefore, the chemical potential affects the binding of the bound states in an opposite way as compared to the momentum space anisotropies. In the semi-QGP region, assuming the potential model is formally unchanged when introducing a chemical potential and the corresponding µ-modifications can be entirely encoded into the Debye mass in exactly the same way as the perturbative case, Eq. (49), then the effective screening mass at non-zero chemical potential reads As a result, one can consider the competition between the two different effects on the binding of a bound state, namely, the anisotropies which lead to a tightly bound quarkonium state and the chemical potential which decreases the binding. In particular, a complete cancellation between the two effects happens when M lm (λ, ξ, µ) = M lm (λ, 0, 0) = Agλ 1 + N f /6. In Fig. 7, the curve on theμ − ξ plane indicates such a complete cancellation which, in the small ξ limit, corresponds toμ/ √ ξ = k(l, m)(N f + 6)/(48N f ). Finally, as mentioned in the introduction, one of the primary motivations for this study was to assess whether it is possible to have an effectively isotropic model potential that reproduces the binding energies of different quarkonium states in a momentum-anisotropic QGP. Our findings can be used to include the effect of momentum anisotropies on in-medium bound state evolution using real-time solution of the Schrödinger equation in a complex screened potential, FIG. 7: The relation betweenμ and ξ. At a given temperature λ, the screening in an anisotropic medium with ξ and µ satisfying this relation is identical to that in an isotropic medium with µ = 0. In this plot, we choose N f = 2. see e.g. [34,35]. However, to do this properly one must prove that the same logic used herein for the real part of the potential can be applied to the imaginary part of the heavy-quark potential. Our preliminary results [36], indicate that the same method can be applied to achieve a numerically reliable isotropic model of the imaginary part of the potential as well. Once this step is complete it will be possible to assess momentum-space anisotropy effects on heavy-quarkonium using isotropic (effectively 1D) simulations.

VI. CONCLUSIONS AND OUTLOOK
In this paper, we introduced a prescription for an isotropic effective Deybe mass that depends on the quantum numbers l and m of heavy-quarkonium state. This mass is obtained by integrating the angle-dependent Debye mass, which emerges when ξ = 0, with spherical harmonic basis functions (3). Tables VII-IX in the appendix demonstrate that when using an isotropic potential with (3) as the effective Debye mass, we can reproduce the energy and binding energies obtained by direct numerical solution of the 3D Schrödinger equation for the same underlying anisotropic potential. Our results demonstrate that, for both small and large anisotropy parameters, one can reproduce the energy to within fractions of an MeV and the binding energy to within a few MeV. As these tables also demonstrate, with this method we are even able to resolve the splitting of the different p-wave states in an anisotropic potential model.
After introducing this method, we discussed why one expects this to be a good approximation even when there is a high degree of momentum anisotropy. We demonstrated that higher-order corrections are under control and that the resulting series converge very quickly, which explains why this prescription does so well in numerically reproducing the full 3D results. Following this, we mentioned some applications of the method introduced herein which include using in real-time solution of the Schödinger equation with a complex in-medium potential. To complete this task, we are now investigating if similar techniques can be applied to the imaginary part of the heavy-quark potential. Preliminary results in this direction are quite promising [36]. small (ξ = 1), moderate (ξ = 10) and large (ξ = 100) anisotropies. Comparing with the results obtained based on the one-dimensional potential model with effective screening masses, the corresponding discrepancies as denoted by δE and δE B are also listed for directly testing our method. For the 3D code, we used a previously developed code called quantumFDTD [10,11,24,[37][38][39].    The exact results of the eigen-energies (E) and binding-energies (EB) for different quarkonium states at various temperatures with ξ = 100. Comparing with the results obtained based on the one-dimensional potential model with effective screening masses, the corresponding discrepancies as denoted by δE and δEB are also listed. All the results are given in the unit of MeV.