Impact of extreme magnetic fields on the QCD topological susceptibility in the vicinity of the crossover region

The topological susceptibility can be formally obtained from the QCD partition function $\mathcal{Z}$ defined through the Euclidean path integral, in the presence of a $\theta$ parameter. The latter is coupled to the topological charge $Q_{\rm top}$ in the Euclidean QCD action $S$ (see e.g. the review Vicari:2008jw ),

where $F_{\mu\nu}^{a}$ is the gluon field strength, $g$ the strong coupling and we also defined the topological charge density $q_{\rm top}(x)$.

Using the functional dependence $\mathcal{Z}(\theta)$, the topological susceptibility is defined via

where $\Omega$ is the space-time volume and $\langle.\rangle$ denotes the expectation value with respect to the partition function $\mathcal{Z}(\theta=0)$. In Eq. (2), we used the translational invariance of the topological charge density correlator as well as the CP-symmetry of the system, ensuring $\langle Q_{\rm top}\rangle=0$. We note moreover that the above definition of $\chi_{\rm top}$ is valid both at zero and nonzero temperature.¹¹1For a discussion on the relationship between the topological susceptibility in Euclidean and Minkowskian space-times, see Ref. Meggiolaro:1998bh .

A homogeneous axion field $a$²²2Only in this section we use $a$ to refer to the axion field. In the rest of the text, it denotes the lattice spacing. couples to gluons in the same manner as the $\theta$ parameter, allowing for the identification $\theta\equiv a/f_{a}$, where $f_{a}$ is the characteristic axion scale. The latter is, in principle, a free parameter of the theory, which can be constrained experimentally Irastorza:2018dyq . Treating $\mathcal{Z}(a)$ as the effective partition function for the axion field as a background field and using Eq. (2), we can see that the susceptibility is proportional to the square of the mass $m_{a}$ of the axion,

This relation is exact and tells us that we can probe the axion mass through considering topological observables in pure QCD. In our study, we will consider QCD with a homogeneous background magnetic field coupled to the fermions. The above definition for the topological susceptibility as well as its relation to the axion mass remains valid while taking into account the relevant, magnetic field-dependent partition function.

In this study we have simulated QCD with $2+1$ flavours of rooted staggered quarks at the physical point with background magnetic fields at finite temperature. The partition function $\mathcal{Z}$ can be written using the Euclidean path integral over the gluon links $U$ as,

where the product in $f$ runs over the up, down and strange quarks and $\beta=6/g^{2}$ denotes the inverse gauge coupling.

The simulations were performed using lattices of geometry $N_{s}^{3}\times N_{t}$, where $N_{s}\,(N_{t})$ denotes the number of points in the spatial (temporal) direction. The physical volume of the lattice is then given by $V=(aN_{s})^{3}$ and its temperature by $T=(aN_{t})^{-1}$, with $a$ being the lattice spacing. We have generated ensembles in lattices of size $24^{3}\times 6,\,24^{3}\times 8,\,28^{3}\times 10$ and $36^{3}\times 12$, varying the inverse coupling in order to scan a range of temperatures around the crossover region, $T=112-212$ MeV, while performing a simultaneous scan for three different magnetic fields, $eB=0,\,0.5,\,0.8$ GeV².

The gluonic action $S_{g}$ is the tree-level Symanzik discretisation, whereas the Dirac operator is discretised using staggered quarks with three times stout-smeared links. The Dirac operator depends both on the gluon links and on the magnetic field, which appears in the product with the charges of the quarks. The latter have been fixed to their physical values $q_{u}/2=-q_{d}=-q_{s}=e/3$, with $e>0$ the elementary electric charge. The masses of the quarks $m_{f}$ have been set to their physical values (assuming degenerate up and down quark masses) using the line of constant physics determined in Ref. Borsanyi:2010cj . The $1/4$ power of the fermion determinant arises due the rooting procedure Sharpe:2006re .

The magnetic field has been considered as a classical background field inside the Dirac operator, hence no dynamical photons were simulated. The electromagnetic potential $A_{\mu}$ was included in a similar fashion as the gluonic one, entering the Dirac operator as $\mathrm{U}(1)$ phases $u_{\mu,f}=\exp{(iaq_{f}A_{\mu})}$ that multiply the $\mathrm{SU}(3)$ gluonic links. The potential was chosen in such a way that it creates a homogeneous magnetic field pointing in the positive $z$ direction. Due to the geometry of our lattices and the periodic boundary conditions for the links, the magnetic field is quantised as $eB=6\pi N_{b}/(aN_{s})^{2}$, with $N_{b}\in\mathbb{Z}$ being the quantum flux number Bali:2011qj . We have to choose $N_{b}$ appropriately for each ensemble in order to maintain the physical field fixed. This implies that there is an uncertainty in the value of the magnetic field in physical units of around 3%.³³3We note that the magnetic field dependence of $\chi_{\rm top}$ could also be determined via a Taylor-expansion around $B=0$. Such Taylor-expansions are possible for magnetic field profiles devoid of the flux quantization condition, see the review Endrodi:2024cqn . The leading-order dependence takes the form of a correlator of $Q_{\rm top}^{2}$ and the magnetic susceptibility Bali:2020bcn , a noisy observable. Since a considerable part of the required $B>0$ ensembles were already generated in Refs. Bali:2011qj ; Bali:2012zg ; Bali:2014kia , we decided to use the direct simulations at $B>0$.

The exact parameters used in the simulations can be found in Appendix A.

The fundamental operator of our study is the topological charge $Q_{\rm top}$. On the lattice, this operator admits several possible discretised forms. We have employed two different lattice definitions for the local operator $q_{\rm top}(x)$ with different scalings to the continuum limit in order to study the systematic errors associated to the choice of discretisation. The simplest discretisation is the one consisting of the smallest closed gluonic loops (plaquettes) and has lattice artifacts of $\mathcal{O}(a^{2})$. Moreover, we also consider an improved version containing larger Wilson loops chosen in such a way that the lattice artifacts of the lowest orders cancel perturbatively. In particular, we use the definition proposed in Bilson-Thompson:2002xlt , which improves the scaling to lattice artifacts of $\mathcal{O}(a^{4})$. However, both discretisations still have lattice artifacts of order $\mathcal{O}(g^{2}a^{2})$. In the following we will refer to the former definition as regular, and the latter as improved.

For a given discretisation of the topological charge density $q_{\rm top}(x)$, the continuum definitions (1) can be simply translated to the lattice. Accordingly, we compute the topological charge and the susceptibility by

where the sum runs over all lattice sites $x$.

Furthermore, in order to calculate the topological charge on a given configuration, we need to remove the ultraviolet fluctuations of the gauge fields. The method chosen for this task is the gradient flow Luscher:2010iy . In particular, we use the Wilson action Wilson:1974sk for the differential equations that define the gauge links at non-zero flow time $\tau_{f}$. For this choice the gradient flow can be also called Wilson flow. Observables constructed from the gauge links obtained for positive $\tau_{f}$ are expected to be renormalised and, thus, finite in the continuum limit Luscher:2013cpa . The gradient flow procedure may also be viewed as an averaging of the fields at each point over a domain with mean-squared radius $\sqrt{8\tau_{f}}$. Hence, only fluctuations with characteristic lengths larger than $\sqrt{8\tau_{f}}$ survive. While this greatly facilitates the determination of the (infrared) topological structure of a configuration, one also has to keep in mind that an excessive amount of flow will destroy the information stored in the gauge fields. At finite temperature, the amount of flow time should not exceed $\tau_{f}^{\rm max}=1/\left(8T^{2}\right)$ Petreczky:2016vrs .

Far from these extremes, the topological charge is expected to be independent of $\tau_{f}$. Hence, an application of the gradient flow to our lattice fields should lead to plateaus as a function of the flow time for the topological charge and susceptibility. We have always made sure that a sufficient amount of flow time has been employed such that those plateaus are reached. In particular, following Ref. Borsanyi:2016ksw , all our observables have been defined at $\tau_{f}=\tau_{f}^{\rm max}$.

To demonstrate that this gradient flow technique is indeed under control, we show the effect of the Wilson flow on our $36^{3}\times 12$ lattices for two different temperatures, above and below the crossover temperature, namely for $T_{1}=150$ MeV and $T_{2}=212$ MeV, with a non-zero background magnetic field $eB=0.5$ GeV², see Figs. 1 and 2. For the integration of the differential equation we have employed a fixed step size of $\Delta\tau_{f}/a^{2}=0.02$.

Once we are close enough to the continuum limit, we are able to obtain almost integer values for the topological charge at all temperatures. We demonstrate this in the left panel of Fig. 1, where we plot the histograms of the topological charge at three different instants of the Wilson flow integration, obtaining almost perfect integers after some finite flow. This can be seen more clearly in the right panel of Fig. 1, where we compare the values of the topological charge for individual configurations before and after the integration. For completeness, we also show how the Wilson flow affects individual configurations, see Fig. 2.

Besides the two different discretisations of the topological charge described above, one can also use alternative definitions that stem from them, namely rounding the value of $Q_{\rm top}$ to the nearest integer or using a different value of $\tau_{f}$, i.e. flowing until 90% of $\tau_{f}^{\rm max}$. This gives eight different definitions, which we studied in our analysis and considered the differences to the improved discretisation using rounding defined at $\tau_{f}=\tau_{f}^{\rm max}$.

Even after employing the Wilson flow for $Q_{\rm top}$, the topological susceptibility suffers from sizeable lattice artifacts that hinder a controlled continuum extrapolation, as shown by previous lattice calculations Borsanyi:2016ksw . Specifically, in our lattice setup $\chi_{\rm top}$ is about an order of magnitude larger than the expected continuum limit at our lowest temperature $T=112\textmd{ MeV}$ and $B=0$ Brandt:2023awt . Most of these lattice artifacts stem from the staggered discretisation of the Dirac operator. Due to the index theorem, on a gauge configuration with topological charge $Q_{\rm top}$, the Dirac operator possesses exactly $|Q_{\rm top}|$ many zero modes (in the continuum)⁴⁴4More precisely, this follows from the index theorem, ensuring that the difference of the numbers of zero modes with positive and negative chirality equals $Q_{\rm top}$ and the fact that in the absence of fine tuning, the Dirac operator only has zero modes with either positive or negative chirality (the so-called vanishing theorem in two dimensions Nielsen:1977aw ; Kiskis:1977vh ). For the lattice equivalent of the index theorem, see also Ref. Hasenfratz:1998ri .. However, the staggered Dirac operator typically does not have exact zero modes. Hence, under the path integral, the weight of each configuration – the determinant of the Dirac operator – is overestimated for topological gauge fields.

To correct for this effect, we introduce the reweighting factor Borsanyi:2016ksw ,

This is nothing but the ratio of the continuum fermion determinant and the staggered determinant in the topological sector, involving the exact zero modes and the corresponding staggered would-be zero modes. In Eq. (6), $m_{f}$ is the mass, $\lambda_{fj}$ is the $j$-th lowest eigenvalue of the staggered Dirac operator for the quark flavour $f$ and $Q_{\rm top}$ is the topological charge, defined via the gluon fields and the Wilson flow, as explained above. The exponent $1/4$ arises due to the rooting procedure and the product in $j$ extends up to the correct number of zero-modes, which for the staggered operator is $4|Q_{\rm top}|$ owing to fermion doubling. Here we also took into account that the staggered eigenvalues come in complex conjugate pairs, and combined the pairs in the second step in Eq. (6).

Therefore, a better estimation of our observable can be obtained via reweighting, i.e. by multiplying each configuration by $W/\langle W\rangle$. This is equivalent to sampling the observable with an improved probability distribution, which is expected to facilitate the continuum extrapolation Borsanyi:2016ksw . Using this method, we can compute the expectation value of a reweighted observable as $\langle\mathcal{O}\rangle_{\mathrm{rw}}\equiv\langle W\mathcal{O}\rangle/\langle W\rangle$. Notice that towards the continuum limit, the staggered would-be zero modes are expected to become exact zero modes Durr:2005ax so that $W$ approaches unity for all configurations. Thus, this kind of reweighting does not change the continuum limit, but it can reduce lattice discretisation errors substantially. The direct effects of the reweighting on the topological susceptibility at finite lattice spacings was reported in Brandt:2023awt .

It is also worth noticing that in Eq. (6) the would-be zero modes are identified with the lowest eigenvalues of the staggered Dirac operator. This identification is clear for configurations at high temperatures, where typical low-lying eigenvalues and would-be zero modes of topological origin are separated in the spectrum. A physical interpretation is more complicated at low temperatures, where topological modes and near-zero modes building up the chiral condensate mix. We will get back to this point below in Sec. 3.2.

We also remark that the larger $|Q_{\rm top}|$, the smaller the reweighting factor, which suppresses the corresponding configuration. Hence, the reweighting in general reduces the topological susceptibility. In turn, for observables that are not directly sensitive to the topological modes, the above reweighting is not expected to have a significant effect. We demonstrate this for the quark condensate $\bar{\psi}\psi_{f}$ in Appendix. B.

The effect of the reweighting is shown in Fig. 3 for our $36^{3}\times 12$ lattice, at a temperature of $T=150$ MeV and with a background magnetic field of $eB=0.8$ GeV². On the scatter plot in the left panel we can observe how for large values of $|Q_{\rm top}|$, the reweighting factor highly suppresses the corresponding configurations. In turn, the contribution of configurations with low topology is effectively enhanced by the reweighting. In the right panel of Fig. 3, we plot the histogram of the topological charge with and without reweighting, where it can be clearly seen how the width of the distribution (the topological susceptibility) is greatly reduced. We note that alternative approaches to reduce staggered lattice artefacts include the use of fermionic definitions for $Q_{\rm top}$ like spectral projectors Athenodorou:2022aay as well as different types of smearing of the link variables, see e.g. Bonati:2014tqa ; Alexandrou:2017hqw .

As commented in the previous section, the low end of the spectrum of the staggered Dirac operator is affected by lattice artifacts. We can distinguish two different regimes for these artifacts, based on the degree of separation between topological and non-topological modes as well as their chirality. The latter is defined⁵⁵5In fact, for our staggered formulation, the correct discretisation of $\gamma_{5}$ is the taste singlet matrix $\Gamma_{5}$, see the definition in e.g. Refs. Durr:2013gp ; Brandt:2024wlw . via the matrix element $\chi_{fj}^{\dagger}\gamma_{5}\chi_{fj}$ for the eigenmode $\chi_{fj}$ corresponding to the eigenvalue $\lambda_{fj}$.

At low temperatures, the would-be zero modes of topological origin get mixed with the low-lying eigenvalues that build up the chiral condensate according to the Banks-Casher relation Banks:1979yr . As the temperature is increased, these low-lying modes disappear and only the topological would-be zero modes remain. We illustrate this picture in Fig. 4. In the left panel of the figure, we see that at low temperatures there is no separation between the low-lying eigenvalues and the topological would-be zero modes, neither in magnitude nor in chirality. However, we see how this situation completely reverses at high temperatures (right panel of Fig. 4) where now the eigenvalues of topological origin are clearly separated from the rest, both in magnitude and in chirality. Accordingly, in this setup one can clearly identify the value of $Q_{\rm top}$ and compare it with the gluonic definition, yielding the same result (as the index theorem dictates). The chirality of the would-be zero modes is not exactly $\pm 1$ due to discretisation effects.

A clear identification of the eigenvalues of topological origin is necessary in order to reweight the fermion determinant, since those modes are the ones that we know become zero in the continuum. As shown in Fig. 4, at high temperatures the identification is clear. However, this becomes problematic at low temperatures, since the would-be zero modes of topological origin may not coincide with the $4|Q_{\rm top}|$ smallest eigenvalues. In order to understand the error introduced by this ambiguity, we have studied the effect of reweighting not the $4|Q_{\rm top}|$ smallest eigenvalues but $4|Q_{\rm top}|$ randomly selected eigenvalues from a window in the spectrum that contains the ${\rm round}(W_{F}\cdot 4|Q_{\rm top}|)$ lowest-lying modes with $W_{F}>1$. We refer to this procedure as multi-range random reweighting and will call $W_{F}$ the window factor.

The reweighted topological susceptibility at finite lattice spacings is found to depend on $W_{F}$ in a non-trivial way. We show this in the left panel of Fig. 5 for all of our lattices at a temperature of $T=112$ MeV with a magnetic field of $eB=0.5$ GeV². This highlights the ambiguity inherent to the reweighting procedure at low temperature and the sensitivity of $\chi_{\rm top}$ to it. However, as we will argue, this effect is largely independent of the magnetic field, so that the reweighting becomes applicable for the ratio of topological susceptibilities at nonzero and zero $B$,

We demonstrate that $R_{\chi}$ (and in particular its continuum limit) is independent of the window size in the right panel of Fig. 5. This tells us that while we cannot distinguish the topological would-be zero modes from the small eigenvalues forming the chiral condensate, the ratios of susceptibilities are insensitive to which particular eigenvalues are chosen. Hence, instead of attempting a continuum extrapolation of the susceptibility on its own, from now on we study the ratio of susceptibilities $R_{\chi}(B,T)$.

We have found that the reweighting procedure is under control for the ratios of susceptibilities and a reliable continuum extrapolation can be performed for them. These results constitute the first non-perturbative determination of the magnetic field dependence of the susceptibility on the lattice. At low temperatures, we can compare the results to chiral perturbation theory estimations, available at next-to-leading order Adhikari:2022vqs . Interestingly, ChPT predicts $R_{\chi}$ to behave in the same way as the similar ratio of chiral condensates Adhikari:2022vqs . We comment on this in Sec. 3.3.1, where we also investigate whether this sum rule is fulfilled beyond the region of applicability of ChPT too. At high temperatures, where the dilute instanton gas approximation is applicable and the temperature scale is much bigger than that of the magnetic field, we expect that the temperature suppression of the susceptibility overcomes any magnetic field dependence, yielding $R_{\chi}\approx 1$. However, the magnitude of the magnetic fields studied in this paper are comparable to the temperatures under consideration. Hence, a non-trivial behaviour is expected even at our highest temperature. We demonstrate that this is indeed the case in Sec. 3.3.2.