Magnetic susceptibility and equation of state of QCD with physical quark masses
Abstract
We determine the free energy of strongly interacting matter as a function of an applied constant and uniform magnetic field. We consider QCD with physical quark masses, discretized on a lattice by stout improved staggered fermions and a tree level improved Symanzik pure gauge action, and explore three different lattice spacings. For magnetic fields of the order of those produced in noncentral heavy ion collisions ( GeV) strongly interacting matter behaves like a medium with a linear response, and is paramagnetic both above and below the deconfinement transition, with a susceptibility which steeply rises in the deconfined phase. We compute the equation of state, showing that the relative increase in the pressure due to the magnetic field gets larger around the transition, and of the order of 10% for GeV.
pacs:
12.38.Aw, 11.15.Ha,12.38.Gc,12.38.MhNow at: ]INFN  Sezione di Pisa, Largo Pontecorvo 3, I56127 Pisa, Italy
Now at: ]School of physics and astronomy, University of Southampton, SO17 1BJ Southampton, United Kindgdom
I Introduction
Strong interactions present some of the most interesting open issues within the Standard Model of particle physics. Quarks and gluons, the elementary colored degrees of freedom of Quantum Chromodynamics (QCD), appear to be confined in ordinary matter, due to a mechanism of color confinement which is not yet fully understood. On the other hand, theoretical arguments cabibbo () and lattice QCD simulations predict the emergence of new phases, in which quarks and gluons are deconfined, when strongly interacting matter is put under extreme conditions, characterized by a high temperature and/or a high baryon density. A strongly interacting plasma of quarks and gluons (QGP), in particular, is the new state of matter which is thought to have filled our Universe in its early stages, and which is hunted for by experiments exploiting ultrarelativistic heavy ion collisions.
We are interested here in the properties of this new “material”, and of strongly interacting matter more in general, when it is placed in the presence of strong magnetic fields. The issue is fundamental, both theoretically and phenomenologically, since in some heavy ion collisions one has the highest magnetic fields ever created in a laboratory hi1 (); hi2 (); hi3 (); hi4 (), reaching up to Tesla (i.e. of GeV), and even larger fields may have been created in the early stages of the Universe vacha (); grarub (). That justifies the recent theoretical investigations on the subject (for comprehensive reviews on various topics, see Ref. lecnotmag ()). In the nonperturbative regime of the theory, this problem can be conveniently approached by lattice QCD simulations buivid (); buivid1 (); buivid2 (); buivid3 (); blum (); reg1 (); reg2 (); levkova (); cffmm (); DEN (); reg0 (); demusa (); thetaeff (); ilgen (); ilgen2 (); reg3 (); reg4 (); bari (). The question regards, in particular, the dependence of the free energy density on a background field , which contains information on the equation of state as a function of and reveals whether the system is diamagnetic or paramagnetic.
The issue is well posed and has been approached by a number of analytic Ioffe:1983ju (); Belyaev:1984ic (); EPS (); Braun:2002en (); Kim:2004hd (); Rohrwild:2007yt (); Goeke:2007nc (); Bergman:2008sg (); agasian08 (); Gorsky:2009ma (); Frasca:2011zn (); F1 (); F2 (); F3 (); endrodihrg (); Nam:2013wja (); Anber:2013tra (); Steinert:2013fza () and recent lattice studies buivid (); reg1 (); reg2 (); cffmm (); levkova (); reg3 (); reg4 (). In a typical lattice setup, the need for magnetic flux quantization introduces a few technical difficulties when one tries to compute the derivatives of the free energy density with respect to . In Ref. cffmm (), we have proposed to overcome such difficulties by determining free energy finite differences, in place of its derivatives, and we have applied this method to the determination of the magnetic susceptibility for QCD, in the standard rooted staggered discretization, with pion masses going down to 200 MeV. Results have provided evidence for a weak magnetic activity in the confined phase and for the emergence of strong paramagnetism in the deconfined, QuarkGluon Plasma phase; such evidence has been confirmed by Refs. levkova (); reg3 (); reg4 ().
The purpose of the present study is to improve on our original determination in several aspects. First of all, we investigate QCD with physical quark masses and three different lattice spacings, in order to check for discretization effects. In this way, we get control over the two major sources of systematic error and provide information which is of direct phenomenological relevance. Second, we will extend the range of temperatures and try to better resolve the region around deconfinement, in order to understand if paramagnetism, or diamagnetism instead, takes place below , and to obtain more information about the high behavior of the susceptibility.
The paper is organized as follows. In Section II we review the method adopted to determine the renormalized free energy density as a function of the background field. In Section III we illustrate the lattice discretization and the numerical setup adopted. In Section IV we present and discuss our results. Finally, in Section V, we draw our conclusions and discuss future perspectives.
Ii The Method
The main focus of our study is on the dependence of the free energy density of QCD on an external constant and uniform magnetic field . Usually, the dependence is given in terms of the derivatives of with respect to , like the magnetization (first derivative) or the susceptibility (second derivative), which are more easily measureable quantities. However, in a lattice setup, taking such derivatives is not trivial, due to field quantization.
Indeed, in order to minimize finite size effects, one usually works on a compact manifold, like a 3D torus. Consistency conditions then require that the magnetic field flux through each closed surface either vanish or be quantized in units of , where is the elementary electric charge of the particles which populate the system. In the case of quarks () moving in a uniform field on a 3D torus one has bound1 (); bound2 (); bound3 (); wiese ()
(1) 
where , are the torus extensions in the directions and must take integer values.
A few approaches have been devised to get around this problem. In Refs. reg2 (); reg4 (), the dependent part of the free energy density has been related to the pressure anisotropy, whose determination requires the knowledge of some perturbative coefficients. In Ref. levkova () the lattice torus has been divided in two identical regions where the magnetic field takes uniform but opposite values, thus ensuring an overall zero magnetic flux without any need for quantization. That allows to take derivatives as usual, even if at the price of introducing interface effects at the boundary between the two regions, which are expected to be negligible in the large volume limit. In Ref. cffmm (), instead, we have developed a method which allows to compute the dependent part of ,
(2) 
where is the partition function of the system and is the spatial volume. As it usually happens when trying to evaluate the ratio of two partition functions umbrella (), the method is computationally demanding, however it is well defined and feasible.
The idea cffmm () is to obtain the finite free energy differences , where and are integers, by integrating the derivative of a suitable extension of the function to real values of ,
(3) 
The operation is well defined as long as the interpolating function is continuous and differentiable, and reduces the problem of determining the ratio of two partition functions to a determination of a standard observable, . It must be stressed that this observable has no direct relation with the magnetization of the original theory, even for integer values of , since it is just the derivative of the interpolating function.
In practice, can correspond to any distribution of magnetic field which interpolates between quantized values. In our realization, the magnetic field will be uniform over the whole lattice apart from a single plaquette, pierced by a sort of Dirac string which brings the excess flux away (see the next section for its explicit form). Of course the final result, , is interpolation independent; that has been also numerically verified in Ref. cffmm ().
Once has been determined, one must take care of getting rid of ultraviolet (UV) divergences, since they contain dependent contributions which do not cancel when taking the difference , and must be properly subtracted. Since we are interested in the magnetic properties of the thermal medium, the most natural prescription is to subtract vacuum (i.e. ) contributions cffmm ():
(4) 
where it is assumed that both terms on the right hand side are computed at the same value of the UV cutoff (lattice spacing). No further divergences are present in Eq. (4), since dependent divergences cannot depend also on , apart from possible finite terms which vanish in the continuum limit (see, e.g., the discussion in Refs. reg0 (); reg2 ()).
The small field behavior of will give access to the magnetic susceptibility of the medium. Indeed, once vacuum contributions have been subtracted, one has the relation
(5) 
where is the magnetization of medium. Assuming that the medium is linear, homogeneous and isotropic, one has , where is the magnetic susceptibility in SI units, so that
(6) 
The field appearing in Eq. (6) is the total field acting on the medium. In our numerical setup, the dynamics of electromagnetic fields will be quenched, so that there is no backreaction from the medium itself, i.e. the magnetization does not change the value of the magnetic field. Therefore, coincides with the applied external field. While that does not introduce any systematic uncertainty in the determination of , which is the physical quantity that we will determine, one should consider that the actual change in the free energy density of a real medium, i.e. capable of producing magnetic fields by itself, will be different. Indeed, if we introduce the auxiliary field , which is generated by external currents only, then
where is the other standard definition of magnetic susceptibility, (). A simple comparison of the two expressions shows that the backreaction of the medium leads to a change of by a factor . However, considering the typical values of that we will show in Section IV (), this correction turns out to be negligible.
In the following it will be convenient to express also in the alternative form
(7) 
which is particularly useful when working in natural units.
Iii Numerical Setup
In this section we discuss our numerical setup for the discretization of quark flavors in the presence of an external magnetic field, with isospin symmetry broken only by the different electric charges. The external electromagnetic field enters the QCD Lagrangian through quark covariant derivatives , where is the abelian gauge four potential and is the electric charge of the quark. On the lattice, covariant derivatives are written in term of the parallel transport ( is the position and the direction) and the introduction of the abelian gauge field amounts to add also an phase: .
The euclidean partition function of the discretized theory in the presence of a magnetic field is expressed as
(8)  
(9)  
(10)  
where is the functional integration over the nonabelian gauge link variables, is the pure gauge action and is the stoutimproved staggered Dirac operator. No integration is performed on the link variables, i.e. the electromagnetic degrees of freedom are quenched.
The action used for the gauge fields () is the tree level improved Symanzik action weisz (); curci (), which involves not only the real part of the trace of the standard square loops () but also that of the rectangles (). The coefficients present in (9) are set to the values and .
In the continuum, the electromagnetic four potential giving rise to a uniform magnetic field can be written, apart from constant terms, in the form and for . The discretization of such a field on a torus is not completely trivial and a possible choice for the phases on the lattice is
(11)  
(12) 
with elsewhere. Here is the Heaviside step function and the term in Eq. (12) is required to guarantee the smoothness of the background field and the gauge invariance of the fermion action (see, e.g., the discussion in Ref. wiese ()).
The particular form of the field given above will be kept also for noninteger values of : that fixes our choice for the free energy density interpolating between quantized values of the magnetic field. This choice is different from other standard discretizations found in the literature (see, e.g., Ref. wiese ()), and in particular from the one used in previous studies DEN (); cffmm (); demusa (); thetaeff (). The difference consists in a simple shift of the Dirac string from one corner to the middle of the lattice. With the present choice, the variation and the magnitude of the phases is minimized, leading to a significant improvement of the signal to noise ratio, as noted in Ref. levkova ().
In the fermionic sector we have adopted the stout link smearing improvement morning () to reduce the effects of the finite lattice spacing and, in particular, to reduce the taste symmetry violations (see Ref. bazavov () for a comparison of the effectiveness of the various approaches). The usual rooting trick is used in Eq. (8) to eliminate the residual degeneracy which is present in the staggered fermion spectrum.
The improved Dirac matrix is built up by means of the two times stout smeared links , which are recursively defined by (see morning ())
(13)  
where the ’s are the original integration link variables, which are used to compute the pure gauge action in Eq. (9). In our simulations we adopted the isotropic smearing parameters . The analyticity of the stout smearing procedure enables to straightforwardly implement the hybrid Monte Carlo update for the improved fermionic action, following Ref. morning ().
Let us explicitly notice that, similarly to Refs. reg0 (); reg1 (); reg2 (), in our simulations the stout smearing is applied only to links, while the external phases are left untouched. See Ref. lecnotmag2 () for a discussion on this point.
24  0.2173(4)  3.55  0.003636  0.1020 

32  0.1535(3)  3.67  0.002270  0.0639 
40  0.1249(3)  3.75  0.001787  0.0503 
We have performed simulations at the physical value of the pion mass, MeV, using the bare parameters , and ( is fixed to its physical value, 28.15) taken from Ref. befjkkrs () and reported in Table 1, which correspond to a line of constant physics at three different values of the lattice spacing .
The number of lattice sites in the spatial direction has been chosen so as to maintain a spatial extent around 5 fm for all values of . At fixed bare parameters, the temperature of the system has been changed by varying the temporal extent of the lattice; in our simulations we explored the temperature range , while our reference symmetric lattices, used to subtract the vacuum contribution, correspond to temperatures below 40 MeV. As we discuss later, in Section IV.1, the fact that the subtraction point is at a low but finite does not introduce any appreciable systematics, due to the rapid convergence to zero of the susceptibility in the low region (exponentially in ).
Finally, the integrand appearing in Eq. (3), , has the following expression:
(14) 
where and are the spatial and temporal extents of the lattice, measured in lattice spacing units. This observable has been measured over a few hundred thermalized trajectories for each parameter set and for each value of , adopting a noisy estimator and averaging over 40 different random vectors for each single measure.
Iv Numerical Results
In Fig. 1 we report an example of the determination of the observable , defined in Eq. (14), for the largest lattice spacing adopted and for two different lattice sizes, and : the former corresponds to MeV, the latter is our reference lattice. The range of explored values of spans the first 5 quanta of magnetic field and for each quantum we have determined on a grid of equally spaced points. Such a grid turns out to be fine enough to allow a reliable integration of (see Ref. cffmm () for a discussion of the related systematic uncertainties).
The integral of over each quantum returns the elementary finite differences . Such quantities are more convenient than the whole difference , since they can be determined independently of each other (one does not need to perform the whole integration from 0 to ) and have therefore independent statistical errors, thus allowing to exploit standard fit procedures for uncorrelated data.
The finite differences corresponding to the data in Fig. 1 are reported in Fig. 2. We also report, after proper rescaling, data obtained on a smaller lattice and at the same value of the lattice spacing, which show the absence of significant finite size effects. Assuming that holds for integer , then
(15) 
Data in Fig. 2 are well reproduced by such a behavior in the whole explored range, and a number of different tests have been performed to check the stability of our fit. In particular, the values obtained for are stable, within errors, if the number of fitted points is changed, and also if a quartic term is added to the free energy, i.e. if a fit of the form is tried, which returns compatible, within errors, with the result obtained by the simple quadratic fit, and compatible with zero. Finally, we have also tried a fit according to a generic power law behavior, , which returns within the precision of .
It is interesting to notice that this implies that strongly interacting matter behaves like a material with a linear response, at least for magnetic fields up to GeV, corresponding to the highest field in the figure. Good linear fits are obtained in similar ranges of for all explored values of and .
The difference of the two slopes, , finally yields the renormalized free energy . The determination of just requires a conversion into physical units for and , according to Eq. (1). The result is
(16) 
in SI units ( and have been reintroduced explicitly). Instead, adopting natural units, one obtains ^{1}^{1}1Actually, and , which are both dimensionless quantities, are related to each other by a simple constant factor, . (see Eq. (7))
(17) 
A similar procedure has been repeated for all combinations of and reported in Table 1. Results are shown in Table 2 and in Fig. 3.
a (fm)  

24  4  0.2173(4)  226(5)  2.648(62)  2.887(68) 
24  6  0.2173(4)  151(3)  0.620(96)  0.67(10) 
24  8  0.2173(4)  113(2)  0.03(10)  0.03(11) 
24  10  0.2173(4)  90(2)  0.02(13)  0.02(14) 
32  4  0.1535(3)  321(6)  4.11(10)  4.48(11) 
32  6  0.1535(3)  214(4)  2.32(12)  2.53(13) 
32  8  0.1535(3)  160(3)  0.86(13)  0.94(14) 
32  10  0.1535(3)  128(2)  0.29(15)  0.32(17) 
32  12  0.1535(3)  107(2)  0.17(15)  0.19(16) 
40  4  0.1249(3)  394(8)  4.62(12)  5.04(13) 
40  6  0.1249(3)  263(5)  2.97(14)  3.24(15) 
40  8  0.1249(3)  197(4)  1.61(14)  1.75(15) 
40  12  0.1249(3)  131(3)  0.27(15)  0.30(16) 
40  16  0.1249(3)  99(2)  0.07(16)  0.07(18) 
iv.1 Discussion
Data displayed in Fig. 3 reveal various interesting features. First of all, one notices that the approach to the continuum limit is very rapid and no significant differences are observed between data computed at different values of the UV cutoff.
Results are in qualitative agreement with those of Ref. cffmm (), obtained by exploiting the same computational method but for with unphysical quark masses, see Fig. 4. Quantitative differences can be explained in part by the presence of the strange quark (which however gives a contribution to the total signal which is not larger than 1/6, see later), and in part by the different pion mass: indeed, an increasing behavior of with decreasing was already observed in Ref. cffmm ().
We thus confirm that strongly interacting matter is paramagnetic, with a magnetic susceptibility which steeply rises crossing the deconfinement transition, which is located around MeV tchot (); tcwup0 (); tcwup1 (); tcwup2 (). However, our data permit to better specify the behavior of the magnetic susceptibility both in the low and in the high region, and in particular around .
In particular the magnetic susceptibility is nonzero and positive, even if relatively smaller, also below , and seems to vanish, within errors, for as low as 100 MeV. It is interesting to notice that data in the low region can be fitted by a simple ansatz , obtaining ^{2}^{2}2In order to check for possible systematic effects related to the choice of the zero temperature subtraction point, which in our case is set to MeV, we have tried to repeat the fit by a function with MeV, but no differences, within numerical precision, have been appreciable, as expected from the fact that the fit returns ., for MeV, MeV, hence in the correct ballpark of the lightest hadrons carrying a nontrivial magnetic moment (starting from the mesons), which are naturally expected to bring the major contribution to the magnetic susceptibility in the hadronic phase. Therefore, this result seems in line with a Hadron Resonance Gas (HRG) model expectation; however, a direct comparison with the HRG prediction of Ref. endrodihrg () is more easily performed in terms of the magnetic contribution to the pressure of the system, and will be presented at the end of this Section.
In Fig. 5 we report a comparison
with other existing lattice determinations
of for QCD.
The results of Ref. reg4 (), which have been obtained
by the same lattice discretization of QCD
but exploiting a different method,
are in quantitative agreement with ours.
Agreement is found also with results
reported in Ref. levkova (), apart from a limited region
below the transition, see Fig. 5.
We do not think this discrepancy to be severe:
it should be reconsidered after continuum extrapolation; moreover,
in the confined region, the finite size effects associated with
the interface in the magnetic field could be more relevant, because
of the larger correlation lengths.
It is interesting to try disentangling the contributions to the magnetic
susceptibility coming from the different flavors.
The observable that we integrate to reconstruct the free
energy is naturally written as the sum of three different
contributions, , see Eq. (14),
and it is therefore straightforward to perform the analysis for the
three different pieces in order to rewrite ^{3}^{3}3
We would like to stress that the separation
into different flavor contributions is not exact,
because of the mixings coming from quark loop contributions.
.
The corresponding quantities are reported separately in Fig. 6,
as a function of temperature. It is maybe more illuminating
to look at the ratios and
, which
are reported in Fig. 7. The former is in nice
agreement, over the whole range of explored temperatures,
with a leading order perturbative expectation based on the squared charge ratio,
. The latter, instead, seems to approach
the same ratio from below in the high limit: this is
expected, since the thermal excitation of
strange degrees of freedom is relatively
suppressed because of the
higher quark mass.
Let us now discuss the behavior of the magnetic susceptibility in the high temperature limit. Our data show an increasing behavior over the whole explored range of temperatures, which however seems to flatten at the highest values of . It is interesting, in this respect, to compare our results with the lowest order perturbative prediction in the regime of asymptotically high temperatures. Based on the results of Ref. EPS (), in the high temperature limit the magnetic susceptibility is given by
(18) 
where and are the quark mass and the quark electric charge in units of and is the number of colors. From the physical point of view, even if the average magnetization of each single particle vanishes in the high temperature limit, because of thermal disorder, the increase in the total density of thermally excited particles and antiparticles compensates that, leading to a logarithmically diverging susceptibility.
It is not reasonable to expect that
our data can be described by the free fermion result, however it is interesting to
notice that values in the
temperature region can be nicely fitted by a logarithmic
behavior, although the coefficient is different from the perturbative one
(see Fig. 8 for a direct comparison of the perturbative
estimate with our results).
Regarding the continuum limit of our results, the density of our data points does not permit a parametrization independent continuum extrapolation and requires instead to fix a definite ansatz for the behavior of . Therefore, we tried to fit our data for the susceptibility in the whole temperature range by a function which reproduces the previously mentioned behaviors in the different regimes:
(19) 
with a continuous and differentiable matching at the temperature , which gives the constraints
(20) 
Although we cannot neglect lattice artefacts (a fit to the whole data would give ) our data are not dense and precise enough to extract the lattice spacing dependence of all the parameters. We explicitly verified that the two sets of fit parameters

with

with
give equivalent results for the continuum extrapolation
(with ).
The value of turns out to be compatible with the transition
temperature, , and the results of the fit
are reported in Fig. 8 and in Appendix A in table format.
One should stress that the errors reported for the continuum extrapolated
values do not take into account the possible systematic error connected to
the particular parametrization chosen for the extrapolation itself.
That could
explain, in particular, the suppression of error bars in the low
temperature region, where the parametrization is exponentially
suppressed, while we do not expect such effect to be significant
in the high temperature region.
Let us discuss now the effect of the magnetic field on the equation of state of the system. The change in the pressure of the system is easily obtained as and is plotted in Fig. 9, as a function of , for two different values of the magnetic field (we make use of our continuum extrapolated determination), normalized by the pressure at (data for the latter quantity have been taken from Ref. presref ()). The different flavor contributions to could be computed as well and, at the quadratic order in , they would be in the same relative ratios shown in Figs. 6 and 7. We notice that the introduction of the magnetic field leads to a relative increase of the pressure which is larger around the phase transition, and in the range 1040% for the typical fields produced in heavy ion collisions at the LHC. In the high regime, instead, the relative increase rapidly approaches zero, as expected, since the pressure at diverges like .
Finally, in Fig. 10, we compare, for a couple of values of , , computed from our continuum extrapolated susceptibility, with the (all orders) HRG prediction for the same quantity, after subtraction of vacuum contributions, extracted from the results reported in Ref. endrodihrg (). A few differences are clearly visible. The HRG model predicts a slightly diamagnetic behavior for low enough , where the contribution from pions, which is indeed diamagnetic, dominates: notice that this behavior of the HRG free energy, which was not underlined by previous literature on the model, regards only the thermal contribution, i.e. after subtracting vacuum contributions to the free energy. For larger temperatures or fields, instead, the contribution of higher spin hadrons becomes overwhelming and the free energy behavior becomes paramagnetic. Lattice results do not confirm this possible weak diamagnetic behavior for small temperatures and fields, even if, for MeV, the possibility is still open, within present statistical uncertainties.
V Conclusions
We have investigated the magnetic properties of strongly interacting matter at thermal equilibrium. The study has been based on a lattice discretization of QCD with physical quark masses: we have considered 3 different values of the lattice spacing and verified the absence of significant discretization effects. We have exploited the method developed in Ref. cffmm (), which is based on a direct computation of the free energy density as a function of a uniform external magnetic background.
We confirm that the strongly interacting medium is paramagnetic, both below and above the deconfinement transition, with a magnetic susceptibility whose order of magnitude is comparable, within the explored range of temperatures, with those of well known strong paramagnetic materials, like liquid oxygen. Moreover, data for the free energy density show that strongly interacting matter behaves like a medium with a linear response, at least for magnetic fields up to GeV, which is the order of magnitude of the typical fields produced in noncentral heavy ion collisions at the LHC; we stress that this behavior is non trivial, since such fields cannot be considered as small, when compared to the typical QCD scale ( GeV).
The magnetic susceptibility is relatively smaller in the confined phase, and compatible with zero within present errors for temperatures MeV. It steeply rises across the deconfinement transition, while in the high regime present data, which go up to MeV, are compatible with a diverging logarithmic behavior, which is predicted by a 1loop calculation of the free energy density.
Notice that the weak diamagnetism, which would be expected at low temperatures and fields based on a pion gas approximation, and which is indeed predicted by HRG computations endrodihrg () (see also Ref. Steinert:2013fza ()), is not confirmed by lattice data, even if present uncertainties still leave room for it for MeV. It would be interesting, in the future, to further investigate this issue.
When compared to the pressure of the standard thermal medium, one sees that the introduction of an external magnetic field leads to a relative pressure increase which rapidly converges to zero in the high phase, while instead it gets larger around the phase transition, where it is in the range 1040% for 0.10.2 GeV, i.e. for the typical fields produced in heavyion collisions. Such contribution becomes rapidly larger and of as the magnetic field increases, so that models for the cosmological QCD phase transition should necessarily take it into account.
We have computed the contributions to the magnetic susceptibility coming from the different flavors: the contribution is approximately of the contribution, as expected on a charge counting basis, over the whole range of explored temperatures; the contribution is slightly smaller than the one of the flavor, but tends to approach it in the high limit.
Future studies, apart from extending the range of explored temperatures and from increasing the precision of data in the low regime, should consider the effects of the inclusion of the quark, whose contribution, especially in the high region, could be comparable to those from the and quarks, since the thermal suppression factor, due to the higher quark mass, should be partially compensated by the electric charge factor. Considering the effects of the inclusion of a baryon chemical potential might be interesting as well.
Finally, we would like to stress that, while present results provide evidence that strongly interacting matter behaves like a medium with a linear response for magnetic fields up to 0.1 GeV, nonlinear corrections could be important for larger fields, which may be of cosmological interest. Our method permits to reconstruct the full dependence of the free energy density on in a straightforward way: while the present study was restricted on purpose to the linear response region, in order to extract the magnetic susceptibility, we plan in the future to extend our analysis to larger fields, so as to determine the nonlinear contributions too.
Acknowledgements.
We thank G. Endrodi and E. Fraga for useful discussions. FS thanks Michael Kruse for useful discussions regarding code optimization on the Blue Gene/Q machine. FN acknowledges financial support from the EU under project Hadron Physics 3 (Grant Agreement n. 283286). Numerical simulations have been performed on the Blue Gene/Q Fermi machine at CINECA, based on the agreement between INFN and CINECA (under INFN project PI12) and on the CSN4 cluster of the Scientific Computing Center at INFNPISA.Appendix A Continuum data for
In Table 3 we report the continuum extrapolated values of , obtained according to the parametrization reported in Eq. (19), for some representative temperatures in the range .
90.0  0.0117(92)  245.0  2.45(10) 
95.0  0.018(13)  250.0  2.53(10) 
100.0  0.028(18)  255.0  2.61(10) 
105.0  0.042(23)  260.0  2.69(10) 
110.0  0.060(30)  265.0  2.77(11) 
115.0  0.084(37)  270.0  2.84(11) 
120.0  0.115(44)  275.0  2.92(11) 
125.0  0.154(50)  280.0  2.99(11) 
130.0  0.203(55)  285.0  3.06(11) 
135.0  0.263(59)  290.0  3.13(12) 
140.0  0.335(60)  295.0  3.20(12) 
145.0  0.420(60)  300.0  3.27(12) 
150.0  0.518(63)  305.0  3.34(12) 
155.0  0.627(68)  310.0  3.40(12) 
160.0  0.743(74)  315.0  3.47(12) 
165.0  0.861(80)  320.0  3.53(13) 
170.0  0.978(84)  325.0  3.60(13) 
175.0  1.094(87)  330.0  3.66(13) 
180.0  1.207(88)  335.0  3.72(13) 
185.0  1.318(90)  340.0  3.78(13) 
190.0  1.426(91)  345.0  3.84(14) 
195.0  1.531(92)  350.0  3.90(14) 
200.0  1.633(92)  355.0  3.95(14) 
205.0  1.733(93)  360.0  4.01(14) 
210.0  1.831(94)  365.0  4.07(14) 
215.0  1.926(96)  370.0  4.12(14) 
220.0  2.019(97)  375.0  4.18(15) 
225.0  2.110(98)  380.0  4.23(15) 
230.0  2.19(10)  385.0  4.28(15) 
235.0  2.28(10)  390.0  4.33(15) 
240.0  2.37(10)  395.0  4.39(15) 
References
 (1) N. Cabibbo and G. Parisi, Phys. Lett. B 59, 67 (1975).
 (2) V. Skokov, A. Y. Illarionov and V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009) [arXiv:0907.1396 [nuclth]].
 (3) V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski and S. A. Voloshin, Phys. Rev. C 83, 054911 (2011) [arXiv:1103.4239 [nuclth]].
 (4) A. Bzdak and V. Skokov, Phys. Lett. B 710, 171 (2012) [arXiv:1111.1949 [hepph]].
 (5) W. T. Deng and X. G. Huang, Phys. Rev. C 85, 044907 (2012) [arXiv:1201.5108 [nuclth]].
 (6) T. Vachaspati, Phys. Lett. B 265, 258 (1991).
 (7) D. Grasso and H. R. Rubinstein, Phys. Rept. 348, 163 (2001) [astroph/0009061].
 (8) D. Kharzeev, K. Landsteiner, A. Schmitt and H. U. Yee, Lect. Notes Phys. 871, 1 (2013).
 (9) P. Cea and L. Cosmai, JHEP 0508, 079 (2005) [arXiv:heplat/0505007]; P. Cea, L. Cosmai and M. D’Elia, JHEP 0712, 097 (2007) [arXiv:0707.1149 [heplat]].
 (10) P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya and M. I. Polikarpov, Nucl. Phys. B 826, 313 (2010) [arXiv:0906.0488 [heplat]].
 (11) P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya and M. I. Polikarpov, Phys. Rev. D 80, 054503 (2009) [arXiv:0907.0494 [heplat]].
 (12) M. Abramczyk, T. Blum, G. Petropoulos and R. Zhou, PoS LAT 2009, 181 (2009) [arXiv:0911.1348 [heplat]].
 (13) P. V. Buividovich, M. N. Chernodub, D. E. Kharzeev, T. Kalaydzhyan, E. V. Luschevskaya and M. I. Polikarpov, Phys. Rev. Lett. 105, 132001 (2010) [arXiv:1003.2180 [heplat]].
 (14) M. D’Elia, S. Mukherjee, F. Sanfilippo, Phys. Rev. D 82, 051501 (2010) [arXiv:1005.5365 [heplat]].
 (15) M. D’Elia and F. Negro, Phys. Rev. D 83, 114028 (2011) [arXiv:1103.2080 [heplat]].
 (16) V. V. Braguta, P. V. Buividovich, M. N. Chernodub, A. Y. .Kotov and M. I. Polikarpov, Phys. Lett. B 718, 667 (2012) [arXiv:1104.3767 [heplat]].
 (17) E. M. Ilgenfritz, M. Kalinowski, M. MullerPreussker, B. Petersson and A. Schreiber, Phys. Rev. D 85, 114504 (2012).
 (18) G. S. Bali, F. Bruckmann, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg, A. Schafer and K. K. Szabo, JHEP 1202, 044 (2012) [arXiv:1111.4956 [heplat]].
 (19) G. S. Bali, F. Bruckmann, M. Constantinou, M. Costa, G. Endrodi, S. D. Katz, H. Panagopoulos and A. Schafer, Phys. Rev. D 86, 094512 (2012) [arXiv:1209.6015 [heplat]].
 (20) G. S. Bali, F. Bruckmann, G. Endrodi, F. Gruber and A. Schaefer, JHEP 1304, 130 (2013) [arXiv:1303.1328 [heplat]].
 (21) C. Bonati, M. D’Elia, M. Mariti, F. Negro and F. Sanfilippo, Phys. Rev. Lett. 111, 182001 (2013) [arXiv:1307.8063 [heplat]].
 (22) L. Levkova and C. DeTar, Phys. Rev. Lett. 112, 012002 (2014) [arXiv:1309.1142 [heplat]].
 (23) M. D’Elia, M. Mariti and F. Negro, Phys. Rev. Lett. 110, 082002 (2013) [arXiv:1209.0722 [heplat]].
 (24) E. M. Ilgenfritz, M. MullerPreussker, B. Petersson and A. Schreiber, arXiv:1310.7876 [heplat].
 (25) G. S. Bali, F. Bruckmann, G. Endrodi and A. Schafer, arXiv:1310.8145 [heplat].
 (26) G. S. Bali, F. Bruckmann, G. Endrodi and A. Schafer, Phys. Rev. Lett. 112, 042301 (2014) [arXiv:1311.2559 [heplat]].
 (27) B. L. Ioffe and A. V. Smilga, Nucl. Phys. B 232, 109 (1984).
 (28) V. M. Belyaev and Y. I. Kogan, Yad. Fiz. 40, 1035 (1984).
 (29) P. Elmfors, D. Persson and B. S. Skagerstam, Phys. Rev. Lett. 71, 480 (1993) [hepth/9305004], Astropart. Phys. 2, 299 (1994) [hepph/9312226].
 (30) V. M. Braun, S. Gottwald, D. Y. Ivanov, A. Schafer and L. Szymanowski, Phys. Rev. Lett. 89, 172001 (2002) [hepph/0206305].
 (31) H. C. Kim, M. Musakhanov and M. Siddikov, Phys. Lett. B 608, 95 (2005) [hepph/0411181].
 (32) J. Rohrwild, JHEP 0709, 073 (2007) [arXiv:0708.1405 [hepph]].
 (33) K. Goeke, H. C. Kim, M. M. Musakhanov and M. Siddikov, Phys. Rev. D 76, 116007 (2007) [arXiv:0708.3526 [hepph]].
 (34) O. Bergman, G. Lifschytz and M. Lippert, JHEP 0805, 007 (2008) [arXiv:0802.3720 [hepth]].
 (35) N. O. Agasian and S. M. Fedorov, Phys. Lett. B 663, 445 (2008) [arXiv:0803.3156 [hepph]].
 (36) A. Gorsky and A. Krikun, Phys. Rev. D 79, 086015 (2009) [arXiv:0902.1832 [hepph]].
 (37) M. Frasca and M. Ruggieri, Phys. Rev. D 83, 094024 (2011) [arXiv:1103.1194 [hepph]].
 (38) E. S. Fraga and L. F. Palhares, Phys. Rev. D 86, 016008 (2012) [arXiv:1201.5881 [hepph]].
 (39) E. S. Fraga, J. Noronha and L. F. Palhares, Phys. Rev. D 87, 114014 (2013) [arXiv:1207.7094 [hepph]].
 (40) J. P. Blaizot, E. S. Fraga and L. F. Palhares, Phys. Lett. B 722, 167 (2013) [arXiv:1211.6412 [hepph]].
 (41) G. Endrodi JHEP 04, 023 (2013) [arXiv:1301.1307].
 (42) S. i. Nam, Phys. Rev. D 87, 116003 (2013) [arXiv:1304.1265 [hepph]].
 (43) M. M. Anber and M. Unsal, arXiv:1309.4394 [hepth].
 (44) T. Steinert and W. Cassing, arXiv:1312.3189 [hepph].
 (45) G. ’t Hooft, Nucl. Phys. B 153, 141 (1979).
 (46) J. Smit and J. C. Vink, Nucl. Phys. B 286, 485 (1987).
 (47) P. H. Damgaard and U. M. Heller, Nucl. Phys. B 309, 625 (1988).
 (48) M. H. AlHashimi and U. J. Wiese, Ann. Phys. 324, 343 (2009) [arXiv:0807.0630 [quantph]].
 (49) G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977).
 (50) P. Weisz, Nucl. Phys. B 212, 1 (1983).
 (51) G. Curci, P. Menotti and G. Paffuti, Phys. Lett. B 130, 205 (1983) [Erratumibid. B 135, 516 (1984)].
 (52) C. Morningstar and M. J. Peardon, Phys. Rev. D 69, 054501 (2004).
 (53) A. Bazavov et al. [HotQCD Collaboration], PoS LATTICE 2010, 169 (2010) [arXiv:1012.1257 [heplat]].
 (54) M. D’Elia, Lect. Notes Phys. 871, 181 (2013) [arXiv:1209.0374 [heplat]].
 (55) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 0906, 088 (2009) [arXiv:0903.4155 [heplat]].
 (56) S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, JHEP 1011, 077 (2010) [arXiv:1007.2580 [heplat]].
 (57) Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B 643, 46 (2006) [heplat/0609068].
 (58) S. Borsanyi et al. [WuppertalBudapest Collaboration], JHEP 1009, 073 (2010) [arXiv:1005.3508 [heplat]].
 (59) A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta and P. Hegde et al., Phys. Rev. D 85, 054503 (2012) [arXiv:1111.1710 [heplat]].
 (60) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabo, arXiv:1309.5258 [heplat].