Deconfinement in the presence of strong magnetic field

We study the impact of a finite magnetic field on the deconfinement phase transition for heavy quarks by computing the fluctuations of the Polyakov loops. It is demonstrated that the explicit Z(3) breaking field increases with the magnetic field, leading to a decrease in the (pseudo) critical temperatures and a shrinking first-order region in the phase diagram. Phenomenological equations that capture the behaviors of the Z(3) breaking field at strong and weak magnetic fields for massive and massless quarks are given. Lastly, we explore the case of dynamical light quarks and demonstrate how an improved constituent quark mass function can enforce the correct magnetic field dependence of the deconfinement temperature in an effective model, as observed in Lattice QCD calculations.

While most effective chiral models (e.g. the NJL model) can capture the effect of magnetic catalysis, they tend to predict the opposite trend in the B-field dependence of the chiral transition temperature [23,24]. To accommodate a decreasing transition temperature with B one usually needs to introduce additional parameter dependence in the effective potentials, e.g. a B-dependent coupling [25], and imposing a specific treatment of the vacuum (and thermal) fluctuations. This is far from ideal as it points to missing interactions (e.g. backreactions, higher order terms, etc.) that are unaccounted for in the original model, demanding a more careful explicit treatment [26][27][28][29][30][31]. 1 1 Even at B = 0, a simple implementation of the PNJL model already leads to a substantially higher transition temperature than LQCD. Implementing a running T d [24,26] may improve the situation, but the curvatures of the potential and various fluctuations remain to be tested.
An important diagnostic test for the correct form of the potential (rather than a parameter change) is to examine higher order fluctuations: as a first step we shall study the Polyakov loop fluctuations in the presence of magnetic field. Ratios of these fluctuations are excellent probes of deconfinement: In a pure gauge theory, they exhibit a jump at the transition temperature [32,33], with well-defined low temperature limits deducible from general theoretical constraints and the Z(3) symmetry. Even in full QCD, they provide a measure for the strength of explicit Z(3) symmetry breaking field induced by the light fermions, and show less renormalization scheme dependence than the Polyakov loop [34,35].
To describe these fluctuations in an effective model, not only the location, but also the curvatures around the minima of the Polyakov loop potential have to be adjusted. The curvatures dictate how reluctant the system is to deviate from the equilibrium position in the presence of external disturbances. Lattice calculations show that the Polyakov loop generally increases with the magnetic field strength [15,17,18,[20][21][22]. Also, the resulting pseudo-critical temperature of deconfinement decreases with the magnetic field [20][21][22], observed for both light and heavy pions [21,22]. On the other hand, the response of the fluctuations to an external magnetic field has not yet been studied, and one of the goals of the current paper is to fill this gap.
In this work we extend the model introduced in Ref. [36] to a finite magnetic field and study its effect on deconfinement in a system with heavy quarks. The study of deconfinement of heavy quarks is interesting in its own right as it relates to other issues such as the modification of the heavy quark potentials in strong fields [37] and color screening [38,39]. In addition, restricting to the heavy quarks allows us to avoid the complications from the chiral transitions. This serves as a foundation for setting up a reliable effective gluon potential to further assess the delicate interplay between light quarks and gluons at finite magnetic field. Comparison between the full quark potential (6) (points), and its linear approximations (8) (lines), versus the real part of the Polyakov loop field for massive (left) and massless (right) quarks. The effect of an external magnetic field on the potential is also shown. For quark masses above 0.8 GeV, the regime explored in this work, the linear approximation provides an adequate description of the full potential.
In our approach the effect of dynamical quarks is modeled by a linear Z(3) breaking term, coupled to the Polyakov loop. This term becomes Landau-quantized when finite magnetic field is present. With such the model setup we explore the Polyakov loop as well as its fluctuations. We shall explore in details how an external magnetic field would enhance the explicit Z(3) breaking, and could eventually change a first order phase transition into a crossover. Moreover, the deconfinement temperature decreases with the magnetic field. We also discuss the phase diagram in the (m l , m s ) plane and find that the external magnetic field shrinks the first order region.
Lastly we explore the case of dynamical light quarks within the PNJL model. The well known problem of an increasing deconfinement transition temperature with B can easily be understood in the current model. We shall also demonstrate, how the use of an improved constituent quark mass function can produce instead the trend of a decreasing deconfinement transition temperature with B.

II. POLYAKOV LOOP POTENTIAL AT FINITE MAGNETIC FIELD
A. robustness of the linear approximation We model the effective Polyakov loop potential as follows: where U G is a pure glue potential and U Q describes the explicit Z(3) symmetry breaking due to dynamical quarks. The Z(3) symmetric pure glue potential U G of choice is from Ref [33]: and A, B, C and D are temperature-dependent model parameters determined from the lattice QCD (LQCD) results on the pressure, the Polyakov loop and its susceptibilities in a pure gauge theory [33]. The U Q potential describes the coupling between quarks and the Polyakov loop. To one-loop order the expression reads [36,40] with In previous studies, we find that the potential U Q can be approximated fairly accurately by the leading term in , i.e.
where the leading constant term U Q ( ,¯ = 0) is dropped as it is irrelevant for computing the Polyakov loop and its fluctuations. The Z(3) breaking strength h Q is given by [36,40] h Q (m, T ) = 12 where K 2 (x) is the modified Bessel function of the second kind. The pure gauge limit is recovered at m → ∞, giving h Q → 0. The opposite limit, m → 0, is also of interest, where we find h Q → 12/π 2 . This means that while the explicit Z(3) breaking field gradually increases as quark masses are decreased, it does not do so indefinitely but saturates at a maximum value. This should be contrasted with the case of chiral symmetry breaking, where the quark mass m serves as the linear breaking field, and can increase without limit. Unfortunately the linear breaking strength for Z(3) is not directly measured in most LQCD studies. A crude estimate [35] based on LQCD results on the ratios of Polyakov loop susceptibilities suggests a larger Z(3) breaking strength than the prediction from the PNJL model. Nevertheless, the issue is far from settled as the analysis is still marred by the unsolved problem of the proper renormalization of the Polyakov loop and its susceptibilities. The study of the quark mass and the mag-netic field dependencies of the Z(3) breaking field may provide some hints for tackling the problem.
The linear approximation made in Eq. (8) is expected to work for heavy quarks. We now examine the efficacy of the scheme for lower (and even vanishing) quark masses. A direct comparison of the full potential U Q and its linear approximation is shown on Fig. 1. Here, the temperature is fixed at T = 0.2 GeV and we show the case of m q = 0.8 GeV (left) and a massless quark (right). We see that the linear approximation provides an excellent description of the full potential for quark masses ≥ 0.8 GeV: the two results are almost indistinguishable. The scheme remains fairly robust even in the worst case scenario of a massless quark, giving only small differences at large values of Polyakov loop x ≥ 0.7.
As a preview we also show in Fig. 1 the influence of a finite magnetic field B on the Polyakov loop potential. The main effect is that the Z(3) breaking strength, corresponding to the slope of the Z(3) breaking potential, is enhanced. We see that for massive quarks even at a moderate qB = 0.4 GeV 2 the linear scheme still provides a very good approximation. In fact for the range of parameters explored in this work, the difference between the full and the linear approximation is negligible. Nevertheless it is not the case for light or massless quarks. For example, we see substantial deviation for massless quark at qB = 0.2 GeV 2 , which is a typical field strength when studying chiral transitions [4].
The versatility of the linear approximation for gaining intuitive understanding shall become obvious in the discussion. In any case, a direct numerical computation of the full potential U Q is rather straightforward and is always an option for quantitative study. In the next sections, a more detailed analysis of the functional dependence on the field strength B and its effects on the fluctuation observables will be explored. In a constant and homogeneous magnetic field background the motion of charged particles undergoes the Landau quantization in the transverse plane. Consequently, the dispersion relation is modified and takes the following form for spin 1/2 particles, where the subsequent Landau levels are quantified by k = 0, 1, 2... and σ = 0, 1 (the spin projection on the B axis). The sum over states is modified accordingly: where the summation runs over Landau levels and the factor |qB|/(2π) accounts for the planar density of each Landau level. Consequently, the explicit Z(3) breaking term becomes a function of the magnetic field, where h B Q (m, T, qB) is obtained by applying the prescription (11) to Eq. (9), The integration over dp z can be performed analytically, which leads to where The main effect of a finite magnetic field on the Polyakov loop is to increase the Z(3) breaking strength.  14), or even the full U Q using the prescription (11), is straightforward and has been explored in previous works. See for example Refs. [4,23] and also Fig. 1. On the other hand, for an intuitive understanding of the magnetic field dependence we shall present both the weak and strong field limits of the linear breaking term, for massive and massless quarks.
We start with the case of massive quarks. Since h B Q is dimensionless, it can be expressed in terms of the dimensionless combinations of m/T and |qB|/T 2 . The weak and strong field limits read: The result is obtained from an asymptotic series expansion of Eq. (14). A detailed derivation, as well as a discussion of higher order terms, are presented in the appendix A. The first term in Eq. (16) coincides with Eq. (9), while the corrections start at the quadratic order, i.e. there is no linear order correction in B in the weak field limit. Due to the coefficient K 0 (m/T ), the response to B is more suppressed for heavier quarks.
When magnetic field becomes strong, the Z(3) breaking field is dominated by the lowest Landau level (LLL), in which while the higher levels are exponentially suppressed. Note, that in the strong field limit the Z(3) breaking field becomes linearly dependent on B. The exact numerical result of Eq. (14) smoothly interpolates between the two limits. See Fig. 2 (right). It is also of interest to study the case of massless quarks. It turns out that the corresponding weak field limit cannot be obtained directly from the m → 0 limit of Eq. (16). Instead we have the following result: For weak fields, beyond the leading quadratic correction, there exist the peculiar odd-power and logarithmic terms. The computational details can be found in the appendix A. Lastly, the strong field limit is dictated by the LLL, and similar to the massive case, the Z(3) breaking field becomes linearly dependent on B.
To summarize the effects of a finite magnetic field to the Z(3) breaking fields: First, it tends to increase the explicit Z(3) breaking strength. Second, the corrections start with quadratic order in qB for weak fields, and are dominated by the LLL at strong fields. The latter gives a linear dependence. Lastly, the increase is more rapid for light than heavy quarks.

III. POLYAKOV LOOP SUSCEPTIBILITIES AND THEIR RATIOS
The Polyakov loop susceptibilities measure the fluctuations of the Polyakov loops. For Z(3) symmetry the order parameter field is complex and one can study fluctuation along the longitudinal and transverse directions. In the language of a potential model, when one focuses on the real sector, these susceptibilities are simply given by the inverse of the curvatures along the real and imaginary directions. The potential U G in Eq. (2) is particularly suited for the current study, as the known susceptibilities at zero explicit breaking, i.e. the SU(3) LQCD results, are reproduced by construction.
To compute the susceptibilities in the current model, we first solve for the expectation value of the Polyakov loop using the gap equation: for φ = (x, y). The susceptibilities are then obtained from the diagonal elements of the inverse of the correlation matrix C: Note, that C should be evaluated using the φ that satisfies the gap equation (19). The susceptibility ratio R T = χ T /χ L can be thus constructed.
The fluctuation of the absolute value of the Polyakov loop, χ A , and the ratio R A = χ A /χ L cannot be obtained simply within a mean-field approach. One way to compute the observable is by setting up a Color Group Integration [35]. Here, we provide an approximate way to determine R A , based on the scaling relation discussed in Ref. [35].
The scaling formula is based on a Gaussian approximation of the potential. In this scheme, R A can be constructed from the mean-field results of 0 , χ L , χ T via: with F given by where K n is the modified Bessel function of the second kind of the n-th order. The scaling variable ξ is computed from 2 : Note, that R A (T, V ) should be computed at a finite volume to be meaningful. Otherwise R A → 1 as V → ∞. In this work, we shall choose V = (6.4 fm) 3 . The numerical results of the Polyakov loop susceptibilities and their ratios will be presented in the next section.

A. fluctuation observables
A key property of a first order phase transition is that it can withstand small external perturbations [41]. For Z(3) symmetry, the transition remains discontinuous for sufficiently heavy quarks, and becomes continuous at the deconfinement critical point. The critical quark mass for a single flavor system in the absence of magnetic field was found to be [36] m 0 = 1.1 GeV.
As shown in Sec II, the presence of an external magnetic field further enhances the explicit Z(3) symmetry breaking. This means that the critical point can be reached at a higher quark mass in the presence of B. To demonstrate this we consider the single flavor case (strange quark) with mass m = 1.4 m 0 , for three values of magnetic field: qB = 0, 3m 2 0 and 6m 2 0 . The observables are presented in Fig. 3

and 4.
We start with the expectation value of the Polyakov loop, shown in Fig. 3 (top left). At B = 0 the Polyakov loop is discontinuous, indicating a first order phase transition. At a sufficiently large magnetic field qB ≈ 3(m 0 ) 2 , the transition becomes continuous, and for an even larger field the transition becomes a crossover.
Similar to the LQCD studies in Refs. [34,42], we can calculate the entropy of a static quark S Q from the Polyakov loop within our model via The peak position of this observable was argued [34] to be a robust way to define the deconfinement transition temperature T D . In particular, T D thus extracted in (2+1)-QCD was found to be lower than that extracted from the inflection point of the Polyakov loop. This trend is also observed in the current model, see Fig. 3 (top right).
In Table I we show the T D 's extracted from the peaks of the observables: χ T , χ L , S Q , ∂ ∂T .
We also check what happens in the case of a larger Z(3) breaking: e.g. for light quarks and/or larger B, and generally find the pattern: T χ T < T S Q T χ L < T inflex. . Differences between characteristic temperatures become substantial in this case.
The change in the order of the phase transition is also evident from the longitudinal susceptibility T 3 χ L , Fig. 4 (left). A key feature of this observable is that it diverges at the critical end point (CEP). This makes the observable χ L unique for defining the (pseudo) critical temperature of deconfinement. Moreover, by studying how the magnitude of the peak changes, the critical value of B can be identified for a given quark mass. On the other hand, the transverse Polyakov loop susceptibility T 3 χ T changes smoothly across the deconfinement critical point, see Fig. 4 (right). In this model it becomes more and more suppressed with increasing B.
In Fig. 5 we show the ratios of the susceptibilities R A and R T . Here, we also show (as points) the results obtained from a Color Group Integration at V = (6.4 fm) 3 . (This volume is sufficiently large to reproduce the meanfield results, e.g. R T .) The observables behave as expected. For R A it interpolates between the two known theoretical limits from 2 − π 2 ≈ 0.43 (26) at low temperatures to 1 at high temperatures. The latter limit is reached more rapidly for a larger breaking strength. We also see that the approximation formula introduced in Sec. III is effective except near the transition points, where we see substantial finite volume effects. For R T , the Z(3) symmetry imposes limit only on the symmetric phase, i.e. 1. The result at high temperatures depends on the model used. Implementing the SU(3) Haar measure is crucial to getting R T 1 at large temperatures, which is also the trend suggested by LQCD studies.
One theoretical motivation for studying and understanding these ratios is to seek out additional observables, other than the standard Polyakov loop, to probe deconfinement. The issue is that the renormalized Polyakov loop extracted by LQCD is multiplied by a scheme dependent (and effectively temperature dependent) factor. It is not clear how to match this quantity to that computed in effective models, and calls into question the physical relevance of the deconfinement features deduced from it, e.g. the transition temperature T D extracted from its inflection point. Similar issue also concerns a direct comparison of the susceptibilities.
It was hoped that by constructing ratios of susceptibil- ities the multiplicative renormalization factor drops out [32]. The task of extracting useful information from these quantities turns out to be more involved than originally thought [33]. In particular, the proper renormalization of the susceptibilities may not be only multiplicative, if true even the ratios would not be scheme independent. As the magnetic field provide an additional handle to study these ratios, the model predictions presented here may guide future LQCD studies at finite magnetic field.

B. phase diagram
To further investigate the effect of a magnetic field on deconfinement we calculated T D by tracking the peak of χ L at various magnetic field strengths. The results are shown in Fig. 6 for the single flavor system, and with quark masses, m s = (0.8, 1, 1.2, 1.4) m 0 .
The deconfinement temperature T D in this model generally decreases when h B Q increases. Thus, when plotted as a function of B, the heavier quarks would require a larger magnetic field strength compared to the light ones, to reach the same T D . This simple observation explains the main features seen in Fig. 6.
Following the arguments in Ref. [36], for effective models of this class, the temperature corresponding to the deconfinement critical endpoint T CEP would remain con- This is verified for the current case of a finite magnetic field, as the B-dependence enters only through the explicit breaking term. Under the same assumptions, the position of deconfinement critical point can be determined from the following condition: where T CEP is given in Eq. (28) and is the critical breaking strength.
We demonstrate two uses of Eq. 29. First, this dictates the critical strength of the magnetic field required to reach the deconfinement CEP for quark masses above m 0 . This is shown by the solid black line on Fig. 7 (left) for the single flavor system, bounded by two limiting lines obtained via the weak and strong field approximations of h B Q in Eq. (16). This agrees with the results obtained by tracking the diverging peaks of χ L , shown as red points in the same figure.
Second, when m s and m l are considered as independent variables, condition (29) determines a critical surface in the (m l , m s ) plane [40]. See also Refs. [43][44][45] for the LQCD determination of this graph. (The continuum extrapolation, however, remains elusive.) For a rough comparison to a recent work [45], the study reported m P S /T CEP = (15.73, 11.15) for N t = (4, 6) in 2-flavor QCD system (NLO). In our 2-flavor model, taking the pseudoscalar mass m P S ≈ 2m crit. , with m crit. ≈ 1.353 GeV and T CEP ≈ 0.261 GeV, the ratio m P S /T CEP ≈ 10.37. This encourages a closer comparative study of effective models and LQCD calculations.
The effect of an external magnetic field on the phase diagram is shown on Fig. 7 (right). Similar to a quark chemical potential [36], an increasing magnetic field tends to shrink the region of the first order phase transition.

V. TOWARDS INCLUDING DYNAMICAL LIGHT QUARKS
It is known that a naive implementation of the PNJL model would give a T D that increases with B. This fact can be easily understood in the current model. 3 We briefly recount the arguments. For potential of the form U = U G − hx, the condition for CEPs reads: The solution fixes the critical values of the Polyakov loop x CEP , hc, and T CEP . If U G is independent of the external fields, the latter two conditions uniquely determine T CEP and x CEP , making T CEP insensitive to the external fields.
Most NJL models can capture the effect of magnetic catalysis, i.e. M Q (T ≈ 0, B) increases with B. This supersedes the enhancing effect of the explicit B dependence on h B Q and leads to a smaller explicit Z(3) breaking in a broad temperature range, including the vicinity of T D . We thus expect the trend of an increasing T D with B for this class of models.
We demonstrate this fact by performing a (2-flavor) PNJL model calculations, based on the NJL model in Ref. [46] and the pure glue potential U G in Eq. (2). The constituent quark mass functions at various values of B, and the Z(3) breaking strength h B Q computed from them, are shown in Fig. 8. As B increases, the diminishing effect of the larger quark mass on h B Q overrides the enhancing effect from the explicit B dependence, suggesting a lower Z(3) breaking. The trend of an increasing T D with B is realized even when the full potential U Q is employed (as in this self-consistent PNJL model calculation), as seen in Fig. 10.
We now explore, under the same mechanism, whether an improved M Q (T, B) would give the correct trend of a decreasing T D with B, as observed in LQCD calculations, e.g. Ref. [22]. To this end, we employ the following parametrization of M Q (T, B): where m = 6 MeV and G ≈ 7 GeV −2 and ψ ψ 0 = −2 × (241.5 MeV) 3 , following the NJL model parameters in Ref. [46]. The essential difference here is, that we employ the lattice result [16] on the ratio of chiral condensates to obtain F: Based on the LQCD results on the quark condensate at zero temperature [16], i.e. F 0 (B) = F(T = 0, B), and the chiral transition temperature at finite B, i.e. T χ (B) [15], we construct a robust parametrization of the function F(T, B): where all quantities are in appropriate units of GeVs. The parametrization is restricted to eB 1 GeV 2 . Note   F(T, B)) from Ref. [16]. that the shape function F 1 can also be well described by a suitably tuned NJL model.
The resulting M Q (T, B)'s are shown in Fig. 9. At low temperatures, the mass functions exhibit a similar increase with B as in the PNJL model case. The key feature, however, is a faster drop of M Q with temperature. This restricts the diminishing effect of a large quark mass, and instead, the enhancing effect from the explicit B-dependence takes over. A direct calculation shows, that indeed h B Q is strengthened by B in the essential temperature range. Employing such h B Q as the explicit Z(3) breaking potential, gives the trend of a decreasing T D with B, see Fig. 10. 4 To examine the dependence on the Polyakov loop potential, we also perform the above analysis for the polynomial potential in Ref. [47]. The result on T D versus B is shown in Fig. 10. When used in a PNJL model, the polynomial potential gives a similar rising trend in T D (B)/T D (0) as U G in Eq. (2). When the improved quark mass function M Q (T, B) (32) is used, the corresponding T D (B)/T D (0) shows the correct decreasing trend for eB < 0.5 GeV, though substantially weaker, and eventually rises again. 5 This clearly demonstrates  (T, B)) and deconfinement. This could help in constraining the missing interactions in effective chiral models. Moreover, a further study could investigate how the explicit and implicit (via M Q (T, B)) B-dependences of h B Q (and the higher order terms) are related to the effects from valence and sea quarks [17]. While the valence quarks always enhance the chiral condensate as B increases, the sea quarks are found to reduce it in temperatures near the chiral transition. Their effects on Z(3) symmetry breaking is left for future research.

VI. CONCLUSIONS
It has been demonstrated that an external magnetic field tends to strengthen the explicit Z(3) breaking. This is a general feature of the one-loop fermionic determinant term. A schematic mean-field calculation shows that the deconfinement phase transition is enhanced, as seen from the lowering of the critical temperature at fixed m or the increase in the critical quark mass. A compact phenomenological formula is derived to capture the effects of a finite magnetic field B on Z(3) breaking: (i) At small B, the correction term ∝ B 2 ; (ii) At large B, the Lowest Landau level dominates and the breaking strength h ∝ B. Extrapolating the results to the case of zero quark masses, we find also terms like B 2 ln B and B 3 .
In addition to quark masses, the magnetic field provides another handle to study Z(3) symmetry breaking. While we demonstrated the inclusion of the higher order terms of the one-loop fermionic potential within meanfield treatments does not lead to substantial changes in the Z(3) breaking strength, we do expect additional modifications from the changes in the Polyakov loop potential, and various beyond mean-field effects, e.g. the spatial dependence of the Polyakov loops and their correlators. We hope that a more detailed study of the fluctuation observables within LQCD can help to distinguish the various effects.
From the perspective of continuum models, the general effects of a finite magnetic field on Z(3) breaking is similar to that of a finite chemical potential. However, for LQCD studies it makes a huge difference: the former would present no sign problem. In particular, we propose the study of the following observables. First, since the logarithm removes the multiplicative renormalization factor from the extracted Polyakov loop. In a continuum model, The second suggestion is to study the magnetic field dependence of the susceptibilities and their ratios, as worked out in this paper. Note, that from the Color Group Integration approach [35] we expect χ A and R A to be explicitly volume dependent. To extract useful information from these quantities, they need to be studied either at a finite volume, or as a function of the scaling variable ξ. Finally, calculations of the correlation function [38,39] (versus distances) and establishing the relation between correlation functions and the susceptibilities will guide the study of the spatial-or momentum-dependence of the Polyakov loop fields. The latter is being explored within the current potential model and will be reported elsewhere. magnetic field, presented in Eq. (16) and Eq. (18).
To facilitate the discussion, we introduce variables x = m/T and dy = 2|qB|/T 2 . The full expression of h Q (x, dy) is given by This expression can be readily computed numerically for arbitrary values of (x, dy). Nevertheless, the strong and weak field limits provide an intuitive way to understand its behavior.

strong magnetic field
The strong magnetic field limit (dy 1) is easy to understand: it is dominated by the Lowest Landau Level (LLL), i.e. the k = 0, σ = 0 part of the sum (A1): Naturally, h LLL depends linearly on dy (the magnetic field). Contributions from the rest of sum (A1) are exponentially suppressed. They appear more prominently as B decreases, more so for the case of small m than for large m. In fact, the correction terms becomes quadratic in the limit of small B, which we turn to next.

weak magnetic field, x = 0 case
We shall employ the Riemann Zeta regularization scheme to derive the expansion of h Q (x, dy) in powers of dy. Noting, that the sum (A1) with σ = 0, compared to σ = 1, differs only by the LLL term, we write The scheme works by directly handling the k-sum in the expansion in powers of dy: (1)] The first term is simply given by the Riemann integral: This recovers the B → 0 result of h B Q in Eq. (9). A key step of the method is to identify the ζ function and represent the various divergent sum by the analytic continuation ζ(−n): These results follow from the general formula where B N 's are the Bernoulli's numbers. It follows, that all sums involving positive even powers vanish as the corresponding Bernoulli's numbers are zero. Note also, that the k = 0 term considered so far does not contribute. It does, however, contribute to the sum Using these results we obtain, up to 6-th order: S(x, dy) = I(x) +c 1 (dy) +c 2 (dy) 2 +c 4 (dy) 4 +c 6 (dy) 6 + . . .
Note, that thec 1 term cancels the LLL term in h B Q in Eq. (A4), giving h B Q (x, dy) = 3 π 2 × (I(x) +c 2 dy 2 +c 4 dy 4 + . . .). (A13) As discussed in the text, the correction term starts from dy 2 order. In Fig. 11 (left) we show the efficacy of the approximation scheme. Note, that Eq. (A13) is an asymptotic expansion, i.e. the series formally diverges, and has a zero radius of convergence. Nevertheless, it can still provide an adequate approximation of the full result using a few terms. The optimal number of terms to keep (N * ) can be inferred by study the ratio r N as a function of N : As the magnitude of the ratio exceeds unity, the accuracy of the sum up to the N-th term, s N , starts to deteriorate. Fig. 12 shows the ratios r N and the corresponding relative errors of s N for typical values of (x, dy). The fact that r N 's are negative means that the series is alternating, also evident from Fig. 11 (left). The value of N * can be extracted where the magnitude of r N exceeds unity. For example, at (x = 2, dy = 10), the series best approximates the exact result when keeping terms up to N * = 4-th order. Including higher order corrections beyond N * would lead to a worse approximation. This also explains the failure of the weak field expansion beyond dy ≈ 12 seen in Fig. 11 (left). Lastly, we note, that N * increases with a heavier quark (large x) or a weaker magnetic field (small dy). In both cases, the accuracy of the scheme improves significantly, as shown in Fig. 12 (right).

weak magnetic field, x = 0 case
The weak field expansion for the x = 0 case cannot be simply obtained from the x → 0 limits of the previous results (c N (x)). For example,c 2 ∝ K 0 (x) → −∞ as x → 0. Nevertheless, both h B Q (x = 0, dy) and S(x = 0, dy) are finite functions of dy. It turns out, that there exists a hidden non-analyticity, hindering the expansion of S(x = 0, dy) in powers of dy.
The first term works out to be I(0) = 4. For the correction terms we see two new complications: First, the existence of terms like (dy) 2 ln(dy), which explains naturally the divergence ofc 2 (x → 0). Second, the contribution from terms such as Again, we observe an explicit cancellation of the linear term, and the correction starts at quadratic order. The effectiveness of the approximation scheme is shown in Fig. 11 (right). Although the accuracy improves with more and more terms, this is of limited use. For large values of dy, the LLL limit h LLL (x = 0, dy) = 3 2π 2 dy, should be used.