Finite-density QCD transition in magnetic field background

Using numerical simulations of lattice QCD with physical quark masses, we reveal the influence of magnetic-field background on chiral and deconfinement crossovers in finite-temperature QCD at low baryonic density. In the absence of thermodynamic singularity, we identify these transitions with inflection points of the approximate order parameters: normalized light-quark condensate and renormalized Polyakov loop, respectively. We show that the quadratic curvature of the chiral transition temperature in the ``temperature--chemical potential'' plane depends rather weakly on the strength of the background magnetic field. At weak magnetic fields, the thermal width of the chiral crossover gets narrower as the density of the baryon matter increases, possibly indicating a proximity to a real thermodynamic phase transition. Remarkably, the curvature of the chiral thermal width flips its sign at $eB_{\mathrm{fl}} \simeq 0.6\,\mathrm{GeV}^2$, so that above the flipping point $B>B_{\mathrm{fl}}$, the chiral width gets wider as the baryon density increases. Approximately at the same strength of magnetic field, the chiral and deconfining crossovers merge together at $T \approx 140\,\mathrm{MeV}$. The phase diagram in the parameter space ``temperature-chemical potential-magnetic field'' is outlined, and single-quark entropy and single-quark magnetization are explored. The curvature of the chiral thermal width allows us to estimate an approximate position of the chiral critical endpoint at zero magnetic field: $(T_c^{\text{CEP}}, \mu_B^{\text{CEP}})= (100(25)\, \text{MeV},\ 800(140)\,\text{MeV})$.

Strongly interacting fundamental particles, quarks and gluons, form a plasma state at sufficiently high temperature. The quark-gluon plasma (QGP), which existed at certain stage of the evolution of the early Universe, may also be created in relativistic heavy-ion collisions. The QGP has been studied at Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory, at the Large Hadron Collider (LHC) at CERN, and will also be subjected to further investigation at Nuclotron Ion Collider fAcility (NICA) at JINR in Dubna, and the Facility for Antiproton and Ion Research (FAIR) in Darmstadt [1].
These experiments offer a unique tool to investigate the QCD phase diagram in a range of increasing baryon densities. A collision of heavy ions creates a QGP fireball which expands, locally thermalizes, cools down, passes through the confining/chiral QCD transition and then (re)hadronizes into final-state colorless states, hadrons. Noncentral collisions also generate a very strong magnetic field which may affect, at least at the early stages, the evolution of the QGP fireball [2]. Despite that the whole process evolves in an out-of-equilibrium regime, certain features of the expanding QGP at zero or sufficiently low baryon density can be determined by its * braguta@itep.ru † maxim.chernodub@idpoisson.fr ‡ andrey.kotov@phystech.edu § amurg@mail.ru ¶ aleksandr.nikolaev@swansea.ac.uk properties in the thermodynamic equilibrium, which are accessible in numerical lattice simulations of QCD.
At vanishing magnetic field and zero baryon density, the equilibrium QCD experiences a broad crossover transition [3] which incorporates the transitions associated with the restoration of the chiral symmetry and the loss of the color confinement in the high-temperature regime.
The crossover has a noncritical character, with both phases being analytically connected. The chiral and deconfining transitions need not to happen precisely at the same point. Moreover, due to the non-singular nature of the crossover, the concrete value of the "pseudo-critical" temperature depends on the operator which is used to define it. The most recent studies indicate that the chiral crossover transition, determined via the inflection point of the light-quark chiral condensate, takes place at T ch c = 156.5(1.5) MeV [4]. The deconfinement transition, identified as the inflection point of the Polyakov loop, appears at substantially higher temperature value, T conf c = 171(3) MeV [3]. Alternatively, one may also use the susceptibilities of these order parameters which would give slightly different crossover transitions even in the thermodynamic limit.
Among many possible options, we define the pseudocritical temperatures of the chiral and deconfining transitions via the inflection points of the light-quark chiral condensate and the Polyakov loop, respectively. These quantities are the order parameters of QCD with quarks of zero masses (the chiral limit of QCD) and with quarks of infinite masses (the pure Yang-Mills theory), where the associated symmetries are not broken explicitly.
Due to the analyticity of the transition, the continuity arguments suggest that the pseudo-critical nature of the transition persists in a low-density region at small values of the baryon chemical potential µ B . Thus, at sufficiently low baryon density, the transition temperature may be expanded over even powers of µ B : where A 2 and A 4 are the first two curvature coefficients of the pseudo-critical transition line. The general form of the polynomial (1) is supported by the analyticity arguments at µ B = 0 along with the invariance of thermodynamic properties of the system under the charge reflection, µ B → −µ B : due to charge conjugation symmetry, the transition temperature of an equilibrium QGP is an even function of the baryon chemical potential µ B . Lattice simulations of the T -µ B phase diagram give the first-principles determination of the transition line (1), which may be confronted with the results of the heavy-ion experiments on the chemical freeze-out line. The freeze-out line corresponds to another curve in the T -µ B plane at which the hadron abundances, that encode the chemical composition of the expanding plasma, get stabilized and thus leave an imprint in the experimentally measured hadronic spectra. It is expected that the chemical freeze-out of the expanding quarkgluon plasma takes place right after the completion of the (re)hadronization process, so that the chemical freezeout temperatures of a generic QGP fireball lies below the pseudo-critical temperature curve (1). The observed momenta of hadrons provide more details on the thermal freeze-out stage that happen at later stages after the chemical freeze-out [5]. The chemical freeze-out temperature may well be described by a polynomial fit similar to the crossover temperature (1) [6].
We study hot strongly interacting matter at low baryonic density subjected to a classical strong magnetic field background. These environmental parameters match the quark-gluon plasma created in the noncentral collisions at the LHC. Due to computational constraints, we do not consider inhomogeneous effects of the high vorticity which is an inevitable feature of plasma created in noncentral collisions with large initial angular momentum [7].
In the first-principles lattice simulations, the effects of the strong magnetic field (B = 0), low baryonic densities (µ = 0) and high temperatures (T ∼ T c ) were studied, so far, in different combinations. At zero magnetic field, the presence of the baryonic matter lowers the pseudo-critical temperature of the QCD crossover transition in the region of low baryon densities. This property is rigidly established in numerical simulations of lattice QCD with imaginary baryonic chemical potential µ I ≡ iµ B [4,8,9] and is also well understood in effective modes of nonperturbative QCD [10][11][12]. A review of recent lattice results may be found in Ref. [13].
At zero baryonic density, the strengthening of the magnetic-field background leads to a smooth decrease of the QCD transition temperature [14]. This phe-nomenon, known as the inverse magnetic catalysis, is not well understood. 1 The strength of the thermodynamic crossover transition was found to increase with the magnetic-field background, possibly indicating the existence of a magnetic-field induced phase transition endpoint at zero baryon density [18]. Notice that the effect of the background magnetic field on transition temperature depends on the masses of the dynamical quarks: at relatively large (unphysical) quark masses, the magnetic catalysis phenomenon is observed: the transition temperature slightly raises with the strength of the magnetic field [19,20] 2 . In this article, we consider QCD with physical quark masses.
Thus, both baryonic density (µ B = 0) and the magnetic field background (B = 0), considered separately, force the temperature of the crossover transition T c to drop. Hence, it would be natural to expect that the combined effect of both these factors, µ B and B, should enhance each other and lead to a much stronger decrease of the crossover temperature.
One of the results of our article is that we confirm the mentioned qualitative expectations. We will also show that the magnetic field affects the magnitude of the leading curvature A 2 = A 2 (B) of the transition temperature T c = T c (µ B , B), Eq. (1). However, we will see that the combined effect of the magnetic field and the baryon density leads to unexpected effects such as strengthening (weakening) of the finite-temperature chiral crossover transition at low (high) magnetic field, with the change of the regime at the magnetic strength of the order of the (vacuum) mass rho-meson squared. An interplay of the wide deconfinement crossover and the narrow chiral crossover is discussed in details. The single-quark magnetization is studied for the first time.
The structure of the paper is as follows. In Section II we discuss particularities of the lattice model and describe technical details of our numerical simulations, which were performed on N t = 6, 8 lattices generated with a Symanzik improved gluons and stout-improved 2+1 flavor staggered fermions at imaginary baryonic chemical potential with subsequent analytical continuation. The properties of the chiral and deconfinement crossovers, uncovered via the chiral condensate of light quarks and the Polyakov loop, are presented, respectively, in Sect. III and Sect. IV. We discuss the pseudocritical temperatures and the thermal widths of both transitions, as well as the effects of the magnetic field and the imaginary chemical potential on these quantities. In Sect. V we use the renormalized Polyakov loop to calculate the single-quark entropy and the single-quark magnetization. The differences between the properties of the magnetization of the bulk quarks and the single-quark magnetization are outlined. The last section is devoted to the discussion of the overall picture of the crossover transitions and to conclusions.

A. Quark densities and chemical potentials
We consider the lattice QCD with three, N f = 2 + 1, quark flavors: two light, up (u) and down (d), quarks and one heavier, strange (s), quark. The total number of quarks N f = d 4 xψ f γ 0 ψ f of the definite flavor f = u, d, s is controlled by the set of the chemical potentials µ f , via the direct coupling in the density part of the action, f µ f N f . The conserved quantities -the baryon number B, the electric charge Q and the strangeness S -are determined by the corresponding chemical potentials µ q with q = B, Q, S, and are related to the quark numbers as follows: Each quark, irrespective of its flavor, carries one-third baryonic charge, and their electric charges are q u = 2/3e and q d = q s = −1/3e, where e = |e| is the elementary charge. The strangeness of the s quark is S = −1, while u and d quarks carry zero strangeness.
Comparing the density part of the action in the basis of the quark numbers and in the basis of the conserved charges, q µ q q = f µ f N f , we find the following relations between all six chemical potentials: In our simulations we take equal potentials for the light quarks and zero chemical potential for the strange quark: With this setup, the chemical potentials for the light quarks, µ, and for the baryon charge, µ B , are related via Eqs. (4), µ B = 3µ. The chemical potential for the electric charge is zero, µ Q = 0. Controversially, the strange chemical potential takes its value from the one of the light quarks, µ S = µ. However, in the absence of a strong electromagnetic background, which would otherwise distinguish between the up and down (strange) quarks due to the difference in their electric charges q f , this choice of the chemical potentials corresponds to near-equal densities of the light-quarks and vanishing strange quark density in the quark-gluon plasma phase. Such quark-gluon plasma should necessarily possess a nonzero (positive) electric charge, but so do the colliding ions. Therefore the choice of the quark content (4) is considered to be a natural one, and it is used in many numerical simulations of the quark-gluon plasma [9,24,25]. In any case, the dependence of the curvature of the phase transition on the chemical potential of the relatively heavy s quark is negligible [26,27], so that we may safely set the chemical potential of the strange quark µ s to zero.

B. Lattice partition function
In our numerical simulations, we partially follow the numerical setup of Ref. [9]. We perform lattice simulations of QCD with N f = 2 + 1 flavors in the presence of purely imaginary quark chemical potentials, µ f = iµ f,I , µ f,I ∈ R, with f = u, d, s, subjected to a strong magnetic field background. We work with the following Euclidean partition function of the discretized theory: where the functional integration is performed over the SU(3) gauge link fields U xµ with the tree-level-improved Symanzik action for the gluon fields [28,29] The lattice coupling β is related to the continuum gauge coupling g in the standard way, β = 6/g 2 . The action (6) is given by the sum over the traces of the flat n × m-sized Wilson lines W n×m x;µν ≡ W n×m x;µν [U ] labelled by the plane vectors µ and ν, and by the starting point x.
The quark degrees of freedom enter the partition function (5) via the product of the determinants of the staggered Dirac operators: constructed from the two-times stout-smeared links U (2) x;ν ≡ U (2) x;ν [U ] following the method of Ref. [30] with the isotropic smearing parameters ρ µν = 0.15 for µ = ν. Here a = a(β) is the lattice spacing. The stout smearing improvement is a standard technique used to ameliorate the systematics related to the effects of finite lattice spacing and reduce taste symmetry violations [31]. Following similar approaches [9, 24-27, 32, 33], we use the rooting procedure in the partition function (5) in order to remove a residual, fourth degeneracy of the lattice Dirac operator (7).
The Dirac operator (7) corresponds to the quarks with the imaginary chemical potential µ f,I subjected to the magnetic field background B. The chemical potential enters the Dirac operator (7) via the additional phases e +iaµ f,I and e −iaµ f,I associated with the temporal links in, respectively, forward and backward directions. The magnetic field appears in the quark operator (7) of the f -th flavor via the composite link field, U f x,µ = u f x,µ ·U (2) x,µ , where U (2) x,µ is the usual (stout-smeared) SU(3) gauge field while the u f x,µ prefactor represents the classical U (1) gauge field corresponding to the uniform magnetic-field background. We consider the classical magnetic background so that the kinetic term of the Abelian field u f x,µ is absent.
In a finite volume with periodic boundary conditions, the total magnetic flux through any lattice plane must be an integer number in units of the elementary magnetic flux [34,35]. For our lattice geometry N 3 s × N t , this property leads to quantization of the strength of the uniform magnetic field B, acting on the quarks of the f -th flavor: Here the integer quantity n ∈ Z counts the number of total magnetic fluxes. Given the fact that the quark electric charges are not the same, one takes the minimal charge, q f ≡ |q d | = e/3, so that the quantization (9) gives a consistent field for all three quarks: For the uniform magnetic field B i = δ i3 B directed along the third axis, the Abelian link field u f x,µ ≡ u f µ (x), acting on the quark of the flavor f , may be chosen in the following explicit form [20]: where x ≡ (x 1 , x 2 , x 3 , x 4 ) is the four-coordinate with the elements running through x ν = 0 . . . N s − 1. The magnetic field B is given by Eq. (9). Due to the periodic structure of the Abelian field (10), the magnetic field cannot be larger then maximal value, determined by the flux number n max = N 2 s /2 , where x gives the greatest integer less than or equal to x. Thus, the nonzero lattice magnetic field B may only be imposed in the range 6π/(eN 2 s a 2 ) B 3π/(ea 2 ), where the strongest value of the field may lead to strong ultraviolet artifacts. To avoid these discretization artifacts, we take n N 2 s /2 in our numerical simulations.

Chiral sector
The chiral condensate ψ ψ is the most straightforward characteristic of the dynamical chiral symmetry breaking in the system of fermions ψ. The condensate vanishes in the phase with unbroken chiral symmetry, ψ → e iγ5ω ψ andψ →ψe iγ5ω , while its deviation from zero signals the violation of the chiral symmetry. The chiral condensate corresponds to an order parameter of the spontaneous chiral symmetry breaking of massless fermions, for which the group of chiral transformations is an exact symmetry group of the classical Lagrangian.
In QCD, the nonzero masses of quarks, m f = 0, break the chiral symmetry explicitly, in all phases. Therefore, the chiral condensate, in a strict mathematical sense, is not an order parameter. However, the condensate of light up and down quarks, with masses well below the characteristic QCD energy scale m u ∼ m d Λ QCD may still serve as an approximate order parameter and thus effectively probe the chiral dynamics.
The chiral condensate of the quark flavor f is given by the partial derivative of the partition function (5) with respect to the quark's mass: where V is the spatial volume of the system. In our N f = 2 + 1 simulations the masses of the light u and d quarks are degenerate, m l ≡ m u = m d . Therefore, it is convenient to introduce the common light quark condensate given by the sum: The chiral condensate of f -th flavor, is evaluated as the trace over the negative power of the Dirac operator (7). Numerically, this calculation is performed with the help of the noisy estimators which comprise O(10) random vectors for each fixed flavor. The finite-temperature renormalization of the lightquark condensate (12) in the presence of the condensate ss of the third, heavier quark s, is implemented following the prescription of Ref. [36]: where m s is the bare mass of the strange quark s. The condensate entering the denominator of the renormalized condensate (14) is computed in the vacuum state, i.e. at zero magnetic field B = 0, zero temperature T = 0, and zero (imaginary) chemical potential µ I = 0.
We took the data for this quantity from (interpolated, when needed) results of Ref. [9]. Other possible renormalization prescriptions may be found in Refs. [24,25].

Gluon sector
The nonperturbative dynamics of the gluon sector gives rise to the confinement of color: the formation of the colorless hadronic states, mesons and baryons, in the lowtemperature QCD. At high temperatures, these states melt, and the system enters the quark-gluon plasma phase with unconfined quarks and gluons. The order parameter of the quark confinement is the Polyakov loop, which may suitably be formulated in the Euclidean QCD as follows: The Polyakov loop operator is averaged over the spatial volume V = N 3 s with the spatial coordinate x. In a purely gluonic Yang-Mills theory, the vacuum expectation value of the Polyakov loop (15) vanishes in the confining, low-temperature phase, and differs from zero in the high-temperature phase that corresponds to the quark-gluon plasma regime. In a purely gluonic theory, the Polyakov loop (15) is an exact order parameter associated with the spontaneous breaking of the global Z 3 center symmetry, P → ZP , where Z = e 2πni/3 , n = 0, 1, 2 are the elements of the center subgroup Z 3 of the SU (3) group. In the presence of light dynamical quarks, the Polyakov loop represents an approximate order parameter of the quark confinement.
For practical reasons of studies of the deconfinement phenomenon, it is convenient to consider the real part of the Polyakov loop:

D. Parameters
We perform numerical simulations at finite temperature around the phase transition using mainly N 3 s ×N t = 24 3 × 6 lattice. In order to estimate the magnitude of the lattice artifacts related to the ultraviolet cutoff effects, we also repeated certain runs on another, 32 3 × 8 lattice with the same ratio N t /N s = 1/4. The comparison of the selected set of results with the ones obtained on the third lattice geometry, 32 3 × 6, gives us an opportunity to estimate the robustness of our data with respect to the finite-volume effects.
The zero-temperature data, used in the renormalization of the condensate (14), were taken from simulations on a 32 4 lattice of Ref. [9]. The physical temperature T = 1/(a(β)N t ) is controlled by the lattice coupling constant β. The lattice spacing varied from a = 0.113 fm at our largest coupling β = 3.7927 till a = 0.253 fm at the lowest coupling β = 3.4949.
The bare (lattice) masses of the quarks, m l and m s , are fine-tuned at each value of the lattice coupling β in order to keep the pion mass at its physical value, m π 135 MeV, and maintain, at the same time, the physical ratio of the quark masses, m s /m l = 28.15. This line of constant physics is well-known phenomenologically from the numerical simulations of Refs. [37][38][39].
We simulated the lattice QCD at the physical point at seven values of the background magnetic field in the interval eB = (0.1 − 1.5) GeV 2 . We took eight points of the imaginary chemical potential of the light quarks, µ I ≡ µ l,I , in the range from a zero value up to µ I /(πT ) = 0.275.
The data for the renormalized chiral condensate may be excellently described by the following function: which has also been used to study the condensate at zero magnetic field in Ref. [9]. The fitting function (17) contains four free parameters: two amplitudes C 0 and C 1 , which describe the scale of the condensate and the degree of its variation over in the crossover region, as well as the pseudo-critical transition temperature T ch c and the width of the crossover δT ch c . All four fitting parameters in Eq. (17) are the functions of the magnetic field B and the imaginary chemical potential µ I .
The numerical data for the condensates and their fits are shown in Fig. 1 for a set of imaginary chemical potentials at smallest nonzero and largest values of the magnetic field. On a qualitative level, the data clearly demonstrate the well-known effect of the inverse magnetic catalysis: the stronger the magnetic field B the smaller the chiral crossover temperature T ch c . They also show that at fixed magnetic field, the increase of the imaginary chemical potential µ I leads, as expected, to increase of the critical crossover temperature.
Before going to the quantitative description of the main results, we estimate the influence of effects of ultraviolet and infrared artifacts of the lattice discretization. At zero magnetic field, these effects were investigated in Refs. [9,26], and we extend the study to the case of the strongest magnetic field, eB = 1.5 GeV 2 , shown in the bottom plot of Fig. 1 for the lowest and largest available imaginary chemical potentials, µ I /(πT ) = 0 and 0.275. The analysis of lattices with different spatial volumes, N s = 24, 32 at fixed temporal extension N t = 6 ensures us that the volume-dependent infrared effects are almost negligible. The inspection of the lattices with N s = 24, 32 and fixed ratio N t /N s = 1/4 demonstrates that while the ultraviolet discretization effects on the condensate are noticeable, the effect of varying lattice spacing on the transition temperature is rather small.
In order to quantify these assertions, we show in Table I the critical temperature T c at both vanishing and largest studied chemical potentials µ I at lattices of all mentioned geometries. The critical temperature, obtained with the fits (17) shown in the bottom plot of Fig. 1, indicate that the variations of the chiral crossover temperature are of the order of 1 MeV, i.e. less than one percent. The lines are the best fits by the function (17). The condensate for the weakest field is shown for all available values of the imaginary chemical potential µI at 24 3 × 6 lattice. The strongest field is represented by the lowest and largest imaginary chemical potentials, µI /(πT ) = 0, 0.275 for 24 3 × 6, 32 3 × 6 and 32 3 × 8 lattices.
We would like to notice that at low (but nonzero) values of the background magnetic fields, there is a particular property of the lattice system which leads to large systematic errors of certain computed quantities. Due to lattice discretization effects, the quantization of the magnetic field (9) limits the number of temperature points at fixed value of magnetic-field strength. Narrowing the study to the crossover region imposes further restrictions thus reducing the quality of the data. We will see below that the data at low magnetic fields has a tendency to possess larger errors at compared to the data at stronger fields.
B. Chiral crossover temperature and its width

General picture
The quantitative analysis of the fits of the chiral condensate gives us the important information how the chiral crossover temperature evolves with increase of the imaginary chemical potential in the magnetic-field background. While the behaviour of the critical temperature is known both at zero chemical potential µ I = 0, Ref. [14] and at zero magnetic field B = 0, Ref. [9], the studies in the full (B, µ I ) plane are performed here for the first time. In addition, we would like to clarify the influence of magnetic field on the thermal crossover width δT ch in the finite-density QCD. This question is important in view of the fact that the role of the magnetic field on the strength of the QCD (phase) transition even at zero chemical potential, µ = 0, has historically been evolving via a set of controversies [14,15,19].
In Fig. 2 we show the spline-interpolated data for the critical temperature of the chiral crossover in the plane of magnetic field B and the squared imaginary chemical potential µ 2 I . One may clearly see that the increase of the imaginary chemical potential, at fixed magnetic field B, leads to the enhancement of the critical temperature for all studied values of B. On the other hand, the strengthening of the magnetic field at fixed imaginary chemical potential µ I gives rise to the decrease of the critical temperature.
The equitemperature curves in Fig. 2 are close to the straight, almost-parallel lines. These properties indicate, respectively, that at small baryon densities (i) the critical temperature of the chiral crossover T ch c at fixed magnetic field B is a quadratic function of the chemical potential µ 2 I ; (ii) the strength of the quadratic dependence does not depend significantly on the strength of the magnetic field. In terms of the baryonic potential, µ 2 B = −(3µ I ) 2 , we conclude that the slope A 2 of the chiral crossover temperature (1) is a positive nonvanishing quantity which moderately depends on the value of magnetic field. The thermal width of the chiral crossover transition (the "chiral thermal width") is encoded in the color of the surface in the same Fig. 2. The chiral width exhibits a weak, but still noticeable, dependence on the imaginary chemical potential. However, the influence of the magnetic field on the chiral thermal width δT ch is much more pronounced: the stronger magnetic field B the narrower transition. This behaviour is well seen in the spline representation of the thermal width in Fig. 3. Interestingly, the magnetic field has a qualitative effect on the behaviour of the chiral width: at weak (strong) magnetic field, the chiral thermal width is an increasing (decreasing) function of the imaginary chemical potential µ I .

Chiral transition temperature and its curvature
At small values of the imaginary chemical potential µ I , the behavior of thermodynamic quantities is necessarily analytic in µ I due to the absence of a thermodynamic singularity in the vicinity the µ I = 0 point. The Taylor series of the observable (real-valued) quantities must therefore run over the even powers of the chemical potential, which makes it possible to use the trivial relation between the imaginary and real baryonic chemical potentials, µ 2 Therefore, the behavior T c = T c (µ B , B) of the critical crossover temperature (1) of the finite-density QCD may be restored from the series of T c (µ I , µ B ) at small imaginary chemical potential µ I : In analogy with the lattice studies with a vanishing magnetic field, we deduce that at B > 0 the curvature A 2 of the critical transition (1) at nonzero baryon density µ B is related to the dimensionless curvature coefficient κ 2 at the imaginary chemical potential µ I in Eq. (18) as: Equations (1), (18) and (19) have rather universal character and can be equally applied to both chiral and deconfining transitions.
In Fig. 4 we show the fits of the critical temperature T c = T c (µ B , µ I ) by the polynomial (18). We fix the magnetic field B and consider the critical temperature as a function of the dimensionless ratio µ I /T . All three fitting parameters T c (B), κ 2 (B) and κ 4 (B) are treated as functions of the magnetic field B. We use both quadratic (with κ 4 ≡ 0) and quartic (with κ 4 being a fit parameter) fits.
The fitting results for the chiral crossover are presented in Fig. 5. We conclude that • The fits allow us to estimate the critical temperature T ch c (B) at zero baryon chemical potential, µ B = 0, subjected to a strong magnetic-field background, Fig. 5(a). The critical temperature decreases with the magnetic field, in an agreement with the inverse magnetic catalysis [14]. In the zero-field limit, our data converge well to the known result T ch c = 156.5(1.5) MeV of Ref. [4], shown by the red square in Fig. 5(a).
• Both for quadratic and quartic fits (18), the quadratic curvature coefficient largely insensitive to the strength of the magnetic field, Fig. 5(b). These fits give qualitatively consistent results, all of which are in agreement with the B = 0 result κ 2 = 0.0132(18) obtained in Ref. [9], and shown by the red square in Fig. 5(b).
• According to Fig. 5(b), the quartic curvature coefficient κ 4 = κ 4 (B) raises with increase of the magnetic field until it reaches the peak around eB fl (0.5 − 0.6) GeV 2 . Eq. (24). At higher magnetic fields, the quartic coefficient κ 4 decreases, and almost vanishes around eB 1.5 GeV 2 . These conclusions have a preliminary character as our numerical results for κ 4 possess rather large statistical errors. Below, we will exclude this coefficient from our analysis, and concentrate on the quadratic truncation of the curvature polynomial (18).
The physical curvature A ch 2 of the chiral crossover temperature (1) for the real-valued chemical potential µ B can be obtained with the help of the analytical continuation (19). The curvature, shown in Fig. 6, seems to exhibit a wide maximum at the magnetic-field strength eB ∼ 0.6 GeV 2 . Unfortunately, the substantial statistical errors of our data do not allow us to determine the presence (and, the position) of this maximum with sufficient certainty. However we will see below that this particular value of the magnetic field marks another interesting effect in the low-density QCD.
To summarize, we observed the effect of the inverse magnetic catalysis both at zero and finite densities. The increasing magnetic field affects the curvature A ch 2 of the chiral crossover transition, making it larger compared to the zero-field value, Fig. 6. We found the presence of the baryonic matter enhances the effect of the inverse magnetic catalysis in a sense that the combined effect of both these factors, µ B and B, leads to a stronger decrease of the crossover temperature.

Chiral thermal width and its curvature
As we have already seen, the thermal width of the chiral crossover transition δT ch c = δT ch c (µ I , B) has a set of interesting features in the parameter space of the magnetic field B and the imaginary chemical potential µ I , as illustrated in Fig. 3. What do these properties mean for the crossover transition in the dense QCD with a realvalued baryonic chemical potential µ B ? In order to answer this question we notice that the thermal chiral width δT ch c -which is, essentially, a difference in temperatures corresponding to opposite sides of the crossover -may be analytically continued to the real chemical potentials, similarly to the critical temperature T ch c itself. To this end, we define the quadratic curvature δκ ch 2 of the chiral thermal width δT ch c as follows: where T ch c (B) ≡ T ch c (µ I = 0, B). Similarity to Eqs. (1), (18) and (19), the thermal width may be analytically continued to the real-valued baryonic potential as follows: is the curvature of the thermal width in the "temperature-baryon chemical potential" plane. In Fig. 7 we demonstrate that the numerical data for the chiral thermal width can be well described by the quadratic function (20). The fits give us the chiral thermal width at zero chemical potential, µ B = 0, shown in Fig. 8(a). The plot suggests that the chiral thermal width is insensitive to the magnetic field until the field reaches the value eB fl ∼ 0.5GeV 2 , and then the width start to decrease slowly.
The effect of the magnetic field background on the curvature of the chiral thermal width is shown in dimensionless, δκ ch 2 , and physical δA ch 2 units in Figs. 8(b) and (c), respectively. It turns out that the curvature δA ch 2 of the chiral thermal width δT ch c may be well approximated by the linear function of the background magnetic field B for both quantities: The best linear fits are shown in Fig. 8(c) by the solid lines. The corresponding best fit parameters are the thermal width at a vanishing magnetic field, δκ   [4,9]. The arrow marks the magnetic flipping point for the width of the chiral crossover (24).
the curvature of the chiral thermal width flips the sign from positive to negative values. We call the value (24) as "the magnetic flipping point". Qualitatively, one can understand the effect of the sign flip of the curvature δA ch 2 as follows. The magnetic flipping point (24) separates two regimes: at weaker magnetic fields, B < B fl , the quadratic curvature is positive, δA 2 (B < B fl ) > 0, and the thermal width of the crossover temperature gets narrower (21) with the rise of the baryon chemical potential. At stronger magnetic fields, B > B fl , the thermal width become wider, δA 2 (B < B fl ) < 0 as the density of the baryonic medium increases.
Numerically, the strength of the magnetic field at the flipping point (24) coincides with the (vacuum) mass of the ρ meson squared, eB fl m 2 ρ 0.601 GeV 2 . At this value of the flipping magnetic field the ρ mesons were proposed to form a superconducting condensate at low enough temperature [40,41]. While this statement is subjected to critical debates [42][43][44], we notice that thermal effects contribute to the ρ meson mass and are likely to destroy the ρ-meson condensate should it be formed at low temperature.
Nevertheless, the closeness of the magnetic flipping point (24) to the mass of the mass scale of the ρ meson suggests that the latter may play a particular role. One could suggest that the mechanism behind the appearance of the magnetic flipping point may be related to the vector meson dominance model [45]. This model proposes that the electromagnetic field interacts with the quark matter via the creation of the quark-anti-quark pairs with the quantum numbers of photons. The lightest such pairs correspond to the neutral rho mesons.
Numerical lattice calculations and effective analytical models suggest that the mass of the neutral meson slowly raises with the strengthening of the magnetic field [42,43,46,47]. Moreover, at the crossover temperature, thermal fluctuations slightly increase the mass of the ρ meson as well [48]. Quantitatively, we expect that the combined temperature and magnetic-field effects at the crossover shift the mass by about 20% from its vacuum value so that the magnetic flipping point (24) is approximately given by the scale of the rho-meson mass.
In order to shed more light on the sign flip of the chiral thermal width, Fig. 8(c), in the next section we study the confining properties of dense QCD in the magnetic field background. The vector dominance hypothesis suggests that the photons interact with the hadronic medium predominantly via the neutral ρ mesons, which are colorless states that do not couple directly to gluons. As we will see below, the sign-flip phenomenon does not exist in the gluonic sector.
We close this section by noticing that our findings on the chiral crossover at zero baryonic density agree well with already known properties of the system. As the strength of the magnetic field increases, the chiral crossover temperature becomes lower, Fig. 5(a) while the transition itself becomes stronger, Fig. 7(a), in agreement with Refs. [14] and [18], respectively.
In addition, our data on the chiral thermal width of the crossover raise an interesting possibility that the parameter plane of the imaginary chemical potential and temperature may contain a thermodynamic phase transition in the limit of large baryonic density at low magnetic field, B < B fl . At stronger magnetic field, B > B fl , the increase of the baryon density leads to the softening of the phase transition.
The chiral pseudo-critical temperature and the thermal width of the chiral crossover, as well the their curvatures in the (µ, T ) plane are summarized in Table II.

Shrinking chiral width and critical chiral endpoint
We would like to finish this section by the following curious observation. As we mentioned above, the width of the chiral transition shrinks in the presence of the baryonic density. Although our numerical simulations are done in the region of relatively low baryon density, we notice that the observation of the shrinking chiral width is consistent with the expectation that at a higher baryon density the crossover turns into a critical endpoint (CEP) of the second order which, at even higher densities, becomes a transition line of the first order.
We may estimate the position of the endpoint as a value of the baryonic chemical potential µ B = µ CEP B at which the width of the phase transition δT ch c (µ B , 0) vanishes. To this end, we use Eq. (21) along with the the results for the chiral thermal width δT ch c (0, 0) and its curvature δA 2 (0) to get for the baryonic density at the CEP: where we neglected the corrections of the order O(µ 4 B ) and higher.
The result (25) is obtained at zero magnetic field. Notice that with the strengthening of magnetic field B, the curvature of the chiral width δA 2 (B) quickly diminishes towards zero, Fig. 8(c), while the width δT ch c (µ B = 0, B) at zero density µ B = 0 drops down less dramatically, Fig. 8(a). Therefore, we may expect that µ CEP B (B) is an increasing function of the magnetic field B. However, at the flipping point (24), our estimation of the CEP (25) formally gives infinite value of the CEP baryonic chemical potential µ CEP B , which shows the limitation of our approach and importance of the higher-order terms, O(µ 4 B ), which were neglected in our simple analysis based on quadratic curvature width (21).
The temperature of the critical endpoint may be obtained with the help of Eq. (1) which takes into account the curvature of the chiral crossover temperature: Curiously, these numbers come quite close to the recent estimation of the location of the critical end point (T CEP , µ CEP B ) = (107, 635) MeV obtained by means of functional renormalization group [49], as well as to other estimations (we refer a reader to Ref. [49] for a detailed review).

A. Renormalized Polyakov loop
The deconfinement (phase) transition is associated with the dynamics of gluons. The corresponding order parameter, in purely gluonic Yang-Mills theory, is the Polyakov loop (15). For the sake of convenience, we study the real part (16) of the Polyakov loop, which is renormalized with the help of the gradient-flow approach following Ref. [50]. The details of the renormalization, and the scheme dependence of the gradient-flow procedure are discussed in the Appendix A. In the vicinity of the crossover, the renormalized Polyakov loop may be well described by the same func-tional behaviour as the chiral condensate (17), where the fitting parameters C 2 and C 3 determine the value of the Polyakov loop at both sides of the crossover region, and T conf c is the pseudo-critical temperature of the deconfinement transition with the deconfining thermal width δT conf c . Some selected fits at lowest (eB = 0.1 GeV 2 ) and highest (eB = 1.5 GeV 2 ) values of magnetic field are shown in Fig. 9 for zero, moderate and highest imaginary chemical potentials, µ I = (0, 0.17, 0.275)πT , respectively.
Similarly to the case of chiral condensate of light quarks, we check the robustness of our results with respect to the volume variations, Fig. 9. In addition of the main results obtained on a 24 3 × 6 lattice, we also show the plots of the renormalized Polyakov loop calculated at a 32 3 × 6 lattice of a higher volume. The visual comparison of the results indicates that the increasing imaginary chemical potential leads to stronger volume dependence at low magnetic field, while at the low density and/or in the strong magnetic field, the sensitivity of the renormalized Polyakov loop to the infrared effects is almost unnoticeable.
The reason for the emergence of these volume effects has a simple systematic origin which is not directly related to the dynamical volume effects. Due to the quantization of magnetic field (9), the number of numerical points, available for the fit (27), is very much limited for weak magnetic fields as compared to stronger magnetic fields. Instead, the magnetic field is varied by the discrete flux variable n = 1, 2, . . . , which needs to be counter-weighted by the variation of the lattice spacing a = a(β) in Eq. (9). The variation of the latter, in turn, affects the temperature T = 1/(aN s ), which quickly goes out of the interesting temperature interval of the deconfining crossover. Therefore we are faced with an artificial limitation of the number of points that could be used in the fit (27), thus bringing a large systematic error to our results. Due to these reasons, we do not discuss below the lattices other than 24 3 × 6 (noticing, at the same time, that the formal low-field B → 0 limit agrees with the known B = 0 results). At larger magnetic fields, the flux variable n may run over larger sets of points and this problem does not exist.

B. Deconfining temperature and its thermal width
In Fig. 10 we show the pseudo-critical temperature of the deconfining crossover as the function of the imaginary chemical potential for the whole set of the available magnetic fields. It turns out that the dependence of the deconfining crossover temperature on imaginary chemical potential can well be fitted by the (quadratically truncated) Taylor series (18) almost at all values of the magnetic field. Notice that the large error bars at the lowest magnetic strength as well as the difference of the results on two lattice sizes N s = 24 and N s = 32 have the systematic origin mentioned above. The temperature T c of the deconfining crossover at zero chemical potential, obtained with the help of the quadratic fits, is shown in Fig. 11(a). Similarly to the chiral crossover temperature, the deconfining crossover temperature is a diminishing function of the magnetic field. This property agrees well with the earlier observation that the gluonic degrees of freedom, as probed by the gluon action, experience the inverse magnetic catalysis similarly to the light quark condensates [51].
According to Fig. 11(a), in the limit of weak magnetic fields, the pseudo-critical temperature of the deconfining transition agrees well with the known B = 0 result, T c = 171(3) MeV, obtained in Ref. [3]. 3 At strong magnetic fields the pseudo-critical line of the deconfining transition, shown in Fig. 11(a), overlaps with the line of the chiral crossover, Fig. 5(a). This fact will be clearer in the last Section, where we discuss the overall phase diagram.
The quadratic curvature of the deconfining transition is shown in Fig. 11(b) as the dimensionless quantity κ conf 2 and in Fig. 11(c) as the physical curvature A conf 2 , calculated via Eq. (19). It seems to have a peak around the magnetic flipping field (24) which is, however, determined with a substantial uncertainty due to large statistical errors. Still, the curvature A conf 2 is a positive quantity so that the pseudo-critical temperature of the deconfining crossover diminishes in the dense QCD matter. This fact means that the presence of the baryon density enhances the "inverse magnetic catalysis" effect for the deconfining phase transition: the presence of matter makes the deconfinement crossover transition happening at lower temperatures.
The thermal width of the deconfining transition may be analytically continued to the real-valued baryonic potential similarly to its chiral counterpart (21). In Fig. 12 we demonstrate that the numerical data for the deconfining thermal width can be well described by the quadratic function (20). In this figure, we dropped a few points with very large error bars which practically do not contribute to the fits while making the figure less readable.
The width of the deconfining crossover is shown in Fig. 13(a). At low magnetic field, the deconfining crossover is very wide, δT conf 60 MeV as compared with the δT ch 11 MeV of the chiral crossover shown in Fig. 8(a). At the magnetic flipping point (24) the width suddenly drops down. At largest studied magnetic field eB = 1.5 GeV 2 the deconfining thermal width δT conf 11 MeV becomes comparable with the chiral thermal width δT ch 5 MeV.
The curvature of the deconfining thermal width is a negatively-valued quantity, as it is shown in dimensionless, Fig. 13(b), and physical, Fig. 13(c), units. The latter has been obtained with the help of Eq. (22) but for the deconfining crossover.
We would like to stress that the presence of the baryonic matter makes the deconfining thermal width wider, thus softening the deconfining transition in the whole studied range of magnetic field.
The deconfining temperature and its thermal width, as well the their curvatures in the (µ, T ) plane are summarized in Table II below. We will discuss the general picture of the chiral and deconfining crossover transitions in the last section.

A. Polyakov loop and thermodynamic potential
The expectation value of the Polyakov loop (15) determines grand-canonical thermodynamic potential of the static quark Ω Q : Free quarks do not exist in the confining phase of QCD. Consequently, at low temperatures, the energy of an individual quark is large and the Polyakov loop is a small quantity. Notice that in the pure Yang-Mills theory the Polyakov loop is an exact order parameter of the quark confinement (the Polyakov loop vanishes in the confinement phase, P = 0) while in QCD the expectation value of the Polyakov loop does not vanish exactly due to the presence of dynamical quarks.
In the deconfining phase the free energy of a single quark is a finite quantity. However, the free energy suffers from ultraviolet divergences due to large perturbative contributions. In order to give a physical meaning to the free energy, it needs to be renormalized. In our paper, we use the gradient flow method to renormalize of the Polyakov loop [50].
In a general thermodynamic system of a volume V , the grand-canonical thermodynamic potential Ω is related to the pressure P as follow, Ω = −P V . The differential of the potential is defined as follows where entropy S, particle number N and magnetization M , determine the response of the grand potential to the variations in temperature T , chemical potential µ and magnetic field B, respectively. These quantities may also be defined for a single static quark introduced into the system by the Polyakov loop operator (15). They have a sense of variation in, respectively, the entropy, the (light) quark number and the magnetization of the overall system in a response of adding one infinitely-heavy, static quark. We use the subscript "Q" for these thermodynamic quantities in order to highlight their single-quark meaning, There is an important feature of our numerical approach: we perform the simulations at a fixed ratio of the imaginary chemical potential µ I to temperature T instead of fixing these quantities separately. Thus for convenience we define the ratios and rewrite the differential of the free energy (29) as Then it is easy to obtain the following relations in terms of (f, T, B) variables: It is worth noticing that the baryon number, determined by a differentiation of the free energy with respect to the baryon chemical potential, is formally an imaginary quantity in our case. However, its physical meaning remains the same, and the baryon number may, in principle, be analytically continued to the domain of the real chemical potential. We do not analyze this quantity in the paper because the accuracy of our numerical data does not allow us to extract the baryon number density unambiguously.

B. Single-quark entropy
Despite the imaginary nature of the chemical potential µ I , the entropy (32) may be determined reliably in the region where the thermodynamic potential is an analytic function of the chemical potential µ. Indeed, Eq. (30) implies the relation f ∂ f ≡ f I ∂ f I which may be used to compute the last term of the entropy of the single quark (32). In this case the entropy of the single quark may be directly expressed via the renormalized Polyakov loop P r : deconfining thermal width where the dimensionless ratio f I is equal to µ I /T .
Although the quark entropy (34) is not sensitive to the Z 3 center symmetry, it is expected to pinpoint a (phase) transition between the low-temperature and high-temperature regions. The entropy of the single heavy quark has a peak -a local maximum as a function of temperature at other parameters fixed -which is close to the pseudo-critical temperature of the chiral crossover [52].
In Fig. 14 we show the single-quark entropy, computed as (34), in the parameter plane of the temperature T and the normalized imaginary chemical potential µ I /(πT ). At each value of the chemical potential the quark entropy has a maximum point which is denoted by a solid blue line. As the imaginary chemical potential increases, the peak is shifted towards higher values of temperature, which is in the qualitative agreement with the picture obtained from our studies of the chiral and deconfinement phase transitions. Unfortunately, with current ensembles we are able to determine the position of the entropy peak only for sufficiently strong magnetic fields, eB > 0.5 GeV 2 . For these large fields, the entropy takes it maximum, in the zero-density limit, at T ∼ 127(10) MeV, in consistency with the position of the common line of chiral and deconfining crossovers.
The calculation of the quark entropy requires an interpolation of the renormalized Polyakov loop to a continuous range of temperatures T and normalized imaginary chemical potentials f I = µ I /T to properly compute derivatives in (34). In this case good enough resolution in terms of discrete (T, µ I /T )-points in the phase transition region is especially important. Unfortunately, we currently have 3 -5 temperature points in the region near T c for all chemical potential and magnetic field values, which turns out to be not enough for the proper estimation of peak in S Q (cubic B-spline with smoothing was employed for interpolation). Moreover, statistics consist of 100 -200 configurations per (T, B, f I ) set, thus relative errors in T c obtained from the maximum of single quark entropy reach 10%. We leave thoughtful study of S Q for future papers.

C. Magnetization
The single-quark magnetization (33) is a real-valued quantity, which may be computed straightforwardly from the renormalized Polyakov loop and then the analytically continued to the real-valued chemical potential. On a first glance the physical meaning of the single-quark magnetization (33) is somewhat obscure since this quantity is associated with the presence of (i) a static, infinite-heavy quark, which (ii) does not possess a spin degree of freedom, and (iii) has zero electric charge.
Due to the latter property, the test quark is not directly coupled to the external magnetic field. Moreover, the immobility of the test quark means that it does not contribute to the Landau diamagnetism, while the absence of the spin, and, consequently, of the magnetic moment, implies the lack of the Pauli paramagnetic contribution. Therefore, one could naively argue that the external test quark would not affect the magnetization properties of the system. On the other hand, the immobile chargeless spinless quark may still affect the electromagnetic properties of the medium since its presence modifiesvia the gluon-mediated interactions -the distribution of the dynamical quarks around it, which, in turn, do couple to the background magnetic field and contribute to the overall magnetization of the system. Therefore, the single-quark magnetization has a meaning of the extent with which the test quark affects the electromagnetically active dense medium of charged quarks and antiquarks.
In Fig. 15 we show the single-quark magnetization (33) computed for three characteristic temperatures at the low (135 MeV), middle (150 MeV), and upper (165 MeV) parts of the crossover transition. The upper row of plots in Fig. 15 corresponds to the actual data obtained for the imaginary chemical potential µ I . In the low-density region, one can perform an analytical continuation of the magnetization data M Im (µ 2 I ) by (i) first expanding the magnetization via the series of the even powers of the imaginary chemical potential µ I ; (ii) and then using the relation, µ 2 B = −(3µ I ) 2 , to get obtain the desired function M (µ B ) = M Im (−(µ B /3) 2 ). At a practical side, we found that the numerical data for M Im (µ I ) at the imaginary chemical potential may be described, at a satisfactory level, by the quartic dependence for all values of magnetic field. Here C i are dimensional fitting parameters. Notice that we make the fits (35) of magnetization at fixed temperatures and magnetic fields so that C i = C i (B, T ). The analytically continued singlequark magnetization, is shown in the lower row of Fig. 15. From these figures one readily observe that the single-quark magnetization is a positive quantity in the whole range of studied parameters (µ I /T, B). In other words, heavy quarks contribute paramagnetically to the overall magnetization of the quark-gluon plasma. Moreover, this paramagnetic contribution is enhanced with the increase of the magnetic field.
Usually, the effect of the heavy quarks on the magnetic polarization of the quark-gluon plasma is ignored because the massive quarks behave as non-relativistic particles for which both the (spin-related) magnetic moment and the (orbital-related) cyclotron frequency are suppressed by the heavy mass. Here we show that (even, infinitely) heavy quarks are magnetically-active constituents of the plasma as they contribute paramagnetically to the overall magnetization.
In order to get a suitable continuous description shown in Figs. 15, we interpolated the data for the single-quark magnetization using the method splines, largely following our approach to of the single-quark entropy S Q . We found that while our data may be used to reliably estimate the inflection point of the Polyakov loop, the presented data for magnetization may contain systematic inaccuracies related to the scarce grid of the data used for the interpolation. This point needs a further investigation.
It is instructive to compare our results on the singlequark magnetization with the behavior of the "bulk" magnetization of the quark-gluon plasma obtained in Ref. [53] in a zero-density limit of QCD. The bulk magnetization of the µ B = 0 quark-gluon plasma is a positive quantity according to [53], thus the zero-density QCD is a paramagnetic medium. The paramagnetic response of the bulk quark-gluon plasma increases in strength both with the increase of temperature and the strengthening of the magnetic field [53]. Our results at µ B = 0 do not show a significant increase in the single-quark magnetization, which may still be consistent with the results of Ref. [53] because our temperature interval (30 MeV) is much shorted compared to that of the quoted reference (about 200 MeV). However, we observe the weakening of the single-quark magnetization with the strengthening of the magnetic field in a sharp contrast with the observed strengthening of the bulk magnetization.
It worth noticing that the single-quark magnetization and bulk magnetization can not be compared with each other directly because these quantities have different physical meanings and they even possess different dimensions: the former quantity corresponds to the energy of a single heavy static quark while the latter number characterizes the energy density of the bulk medium. Nevertheless, both quantities characterise the magnetic properties of the strongly interacting medium subjected to an intense magnetic field background.

VI. OVERALL PICTURE AND CONCLUSIONS
In our work, we studied the influence of the strong magnetic field on the chiral and deconfinement transitions in finite-temperature QCD at a low baryonic chemical potential. In the low-density QCD with real (physical) masses of u, d, and s quarks, these transitions are not accompanied by any thermodynamic singularities in the parameter space of the theory. Instead, the theory experiences a smooth broad crossover from the cold chirallybroken hadronic medium to the hot chirally-symmetric plasma of deconfined quarks and gluons.
In the absence of a real phase transition, the positions of the chiral and deconfining crossovers are not well defined; they depend on a particular form of the operator employed to probe them. In our paper, we identify the location of the chiral crossover as the inflection point of the expectation value of the chiral condensate of light quarks, which is an exact order parameter for the chirally broken phase in QCD with massless quarks.
We reveal the location of the deconfining crossover via the inflection point of the expectation value of the gradient-flow-renormalized Polyakov loop, which is the order parameter for the deconfinement phase transition in a pure Yang-Mills theory (QCD with infinitely massive quarks).
In addition to the positions of the chiral and deconfining crossover lines in the parameter space, we determined the thermal width of each of these crossovers. In the absence of thermodynamically singular behaviour, the thermal width may serve as a quantitative characteristic of the strength of the crossover transition. The thermal width δT , formally defined via the fitting functions (17) and (27) -can be understood as a temperature range over which the corresponding order parameter reaches its values at both sides of the crossover. Similarly to the positions of the crossover lines, their thermal widths are prescription-dependent quantities which may depend on the operator used to identify them.
We performed the calculations on N t = 6, 8 lattices generated with Symanzik improved gluons, and stoutimproved 2+1 flavor staggered fermions at physical quark masses and imaginary baryonic chemical potential. We used the analytical continuation from purely imaginary to real-valued baryon chemical potential.
Below, we summarize all effects of the magnetic-field background on the chiral and confining crossover transitions at low baryonic densities and finite temperature.
The chiral crossover: 1. The effect of the inverse magnetic catalysis extends to the region of low baryon densities: the chiral crossover temperature drops down as the background magnetic field strengthens. Moreover, the presence of the baryonic matter, the effect of the inverse magnetic catalysis becomes slightly stronger.   (27), respectively. We show the pseudo-critical temperatures Tc, the widths δTc, as well as the dimensionless quadratic curvatures of the crossover temperature κ2 and its width δκ2, determined, correspondingly via the quadratically truncated fits (18) and (20). The curvatures in the physical units, A2 and δA2, are found via the fits (19) and (22). The marks denote the data taken from other sources: (a) Ref. [4], (b) Ref. [3] and (c) Ref. [9]. The data point (d) is derived from (c) via Eq. (19). Note that for our data we present only the statistical errors, while the points (a)-(d) include also systematic uncertainties coming from an extrapolation to the continuum limit.  (38) changes its sign at the magnetic flipping point, eB 0.6 GeV 2 , as shown in Fig. 8(c). Thus, the presence of the baryon matter makes the chiral crossover transition narrower (wider) in the magnetic-field background with B < B fl (B > B fl ).
8. The curvature of the deconfining transition is a positive-valued quantity with a peak around the critical value of magnetic field eB fl , Fig. 11(c). Therefore, the presence of the baryonic matter lowers the deconfining temperature at finite magnetic field, thus, effectively, enhancing the inverse magnetic catalysis of the deconfining crossover. The maximum enhancement happens around the magnetic flipping field B B fl , Eq. (24).
9. The deconfining crossover is generally a much wider transition as compared to the chiral crossover. This fact follows from comparison of their widths, Fig.  8(a) and Fig. 13(a), respectively. However, the deconfining thermal width decreases very rapidly with magnetic field: it shrinks at least five times, from δT conf c 60 MeV at a vanishing field to δT conf c 11 MeV at the strongest studied field, eB = 1.5 GeV 2 .
10. The curvature of the thermal width of the confining crossover is a negative quantity in the whole studied range of magnetic fields, Fig. 13(c). It means that the presence of baryonic matter tends to weaken the deconfining crossover in the studied range of B.
The overall picture of the crossover region: The pseudo-critical temperatures and the thermal widths of the chiral and deconfining crossovers, as well as the their curvatures in the (µ, T ) plane, are summarized in Table II. We illustrate the overall picture of the crossover region in the (B, T ) plane in Fig. 16. We show the pseudocritical temperatures T c of the chiral and deconfining transitions, as well as their thermal widths δT c as functions of magnetic field B for three different values of the baryonic chemical potential: µ B = 0, 250, 500 MeV. We used the quadratic analytical continuation for the chiral crossover temperature (37), its width (38), and the same quantities for the deconfining crossover transition. The value of the largest chemical potential, µ B = 500 MeV, is specially chosen for illustrative purposes in order to highlight the qualitative effects of the baryonic matter on the phase transition. At this relatively high baryon density, the presence of the quartic term in the Taylor expansions over the chemical potential may affect, quantitatively, both the transition lines and their widths.
The three plots in Fig. 16 capture all basic properties of the crossover transition region: 11. At a vanishing magnetic field and zero baryonic density, the deconfining transition is a wide crossover with the thermal width δT conf 60 MeV which takes place at T conf 170 MeV. The chiral transition is much narrower crossover, δT ch 11 MeV, that takes place at somewhat lower temperature, T ch 156 MeV.
12. As the magnetic field strengthens, the transitions temperatures of the confining and chiral crossovers become lower (the inverse-magnetic catalysis phenomenon). Both crossover transitions become narrower and, therefore, stronger.
13. At the "tri-pseudo-critical" point (eB * , T * ) (0.5 GeV 2 , 140 MeV) these transitions merge together and overlap at higher magnetic fields, with different widths for chiral and deconfining crossovers. The tri-critical point appears at the magnetic flipping field, B * B fl . Since the deconfining crossover is very wide (δT conf > |T conf c −T ch c | in the whole studied region), the merging point has a rather academic significance.
14. The presence of the baryonic matter enhances the inverse-magnetic-catalysis effect for both crossovers: the pseudo-critical temperatures drop as the baryon chemical potential increases in the whole studied region of magnetic field. Both chiral and deconfining curvatures are found to be (generally) increasing in the presence of the magnetic field.
15. The presence of the baryonic matter always widens the deconfining crossover. 16. The effect of the baryon density on the chiral crossover is two-fold: the matter makes chiral transition narrower (wider) at lower (higher) fields compared to the magnetic flipping field (24) B fl ≈ B * , where the both crossovers merge.
17. The behaviour of the chiral thermal width and its curvature give a simple estimation (26) of the critical endpoint in the T − µ plane in, surprisingly, reasonable range of parameters (26): ) (100, 800) MeV.
In addition, we have studied the single-quark entropy and the single-quark magnetization. The maximum of the single-quark entropy (32) corresponds very well to the common chiral-deconfining transition line at larger magnetic fields, B B * , Fig. 14. Since the deconfining crossover is very wide, thus the peak of S Q is broad and it is difficult to pinpoint the maximum of the entropy with acceptable accuracy at lower magnetic fields with our current statistics.
The single-quark magnetization (33) exhibits a variety of nontrivial features, Fig. 15. First of all, the magnetization of the heavy quarks turns out to be nonzero. This fact reveals a surprising property of the system because for non-relativistic heavy particles both the magnetic moment and the cyclotron frequency are suppressed by the large mass, implying -naively -that these particles do not contribute to magnetic properties of the plasma. We argue that the influence of the heavy quarks on the magnetization goes indirectly. The heavy quarks affect locally the dynamics of the light quarks, while the latter quarks, being magnetically active, contribute to the excess of the overall magnetization.
Given the scarce number of points in the direction of magnetic field, the following features of the single-quark magnetization -made in the vicinity of the crossover transition at T = (135 − 165) MeV -are largely of a qualitative nature: (i) The single quark has the paramagnetic response to the external magnetic field (i.e., the single-quark magnetization is a positive quantity at all studied fields).
(ii) As the strength of the magnetic field increases, the magnetization drops down in the low-density crossover region.
(iii) At zero magnetic field and low baryonic densities, the magnetization, as the function of the realvalued baryon potential, increases (decreases) in the hadronic (quark-gluon plasma) regions of the crossover.
It is important to stress that the exact positions of the crossover transitions and their thermal widths are prescription-dependent quantities. Their precise values depend not only on the operators used to reveal them but also on the (re)normalization of these operators. However, the analysis indicates that our results match well with other available data at the corners of the explored parameter space, thus providing us with additional support for the validity of the presented picture in the whole explored region. In this paper we renormalize the expectation value of the Polyakov loop using the gradient flow procedure [54,55]. This procedure removes perturbative ultraviolet content of gauge fields. We refer an interested reader to Refs. [50,54] for a comprehensive account of the renormalization with the help of gradient flow.
The approach postulates the evolution of the gauge field configurations in the "Wilson flow" 4 space defined by the differential equation of a diffusion type: The flow time τ controls the degree of smoothing of the initial gauge configuration (A2) in the space of gaugefield configurations, guided by the gauge action S YM and the bare gauge coupling g 0 . The functional derivative ∂ x, µ in Eq. (A1) acts in the Euclidean coordinate space and in the color space [54]. The dot over V in Eq (A1) denotes a partial derivative with respect to the flow time τ . Due to the diffusive character of the evolution equation (A1), the flow smears gluon configurations at the length scale The operators built from the flow-evolved variables V x, µ do not require an additional renormalization [56] after extrapolation to τ → 0 limit. In the case of Polyakov loop for small τ , the evolved operator is equivalent, up to a multiplicative factor, to the renormalized original operator at τ → 0 [55,56]. The gradient flow procedure has an intrinsic ambiguity related to the choice of the "optimal" flow time at which the smoothing procedure should stop. On the physical grounds, the optimal τ is naturally constrained within the ultraviolet and infrared limits, a √ 8τ Λ QCD . However, this interval is too broad to fix the renormalized free energy (28) Ω Q unambiguously. Since the Polyakov loop is renormalized multiplicatively, the uncertainty in the renormalization scale leads to an additive ambiguity in the renormalized free energy determined at two scales f andf , related to the corresponding optimal flow times via Eq. (A3). In order to probe the dependence of the free energy on the choice of the optional flow time, we compared the shift (A4) for two different values of the renormalization scale,f = 0.54 fm and f = 0.80 fm, for two values of magnetic field, eB = 0.5 GeV 2 and eB = 1.5 GeV 2 . The results are shown in Figs. 17(a) and (b).
It appears that the energy shift (A4) does not depend, within the error bars, on the imaginary chemical potential, i.e. chemical potential does not affect the renormalization. For a moderate magnetic field eB = 0.5 GeV 2 , the energy shift ∆Ω Q is almost a temperatureindependent quantity, with a slight systematic dropalbeit within the large error bars -at the colder side of the pseudo-critical crossover temperature T 140 MeV, see Fig. 17(a). On the contrary, according to Fig. 17(b), the drop in the shift at the low-temperature side of the crossover region is clearly visible at the stronger field, eB = 1.5 GeV 2 . Figure 17 suggests that the energy shift ∆Ω Q for different renormalization scales f andf is affected by the strong magnetic field background. Let's estimate if this scheme dependence may influence the determination of the deconfining crossover temperature from the inflection point of the Polyakov loop. The energy shift δΩ Q differs from ∼ 50 MeV at the cold (T − ∼ 100 MeV) side of the crossover to ∼ 55 MeV at the hot (T + ∼ 200 MeV) side. Thus exp −∆Ω Q T −1 + − T −1 − ∼ 1.3, i.e. the multiplicative bias in renormalization is about 30% in both ends. But the change in the magnitude of the Polyakov loop, induced by the deconfinement phenomenon (Fig. 9), amounts to the factor of 10, which is about 30 times bigger compared to the mentioned systematics of the scheme. Thus we expect that this effect of the renormal-ization scheme dependence may be safely neglected. We also expect that the single-quark entropy and the singlequark magnetization are not affected by the described gradient flow renormalization systematics.