Meson masses in electromagnetic fields with Wilson fermions

We determine the light meson spectrum in QCD in the presence of background magnetic fields using quenched Wilson fermions. Our continuum extrapolated results indicate a monotonous reduction of the connected neutral pion mass as the magnetic field grows. The vector meson mass is found to remain nonzero, a finding relevant for the conjectured $\rho$-meson condensation at strong magnetic fields. The continuum extrapolation was facilitated by adding a novel magnetic field-dependent improvement term to the additive quark mass renormalization. Without this term, sizable lattice artifacts that would deceptively indicate an unphysical rise of the connected neutral pion mass for strong magnetic fields are present. We also investigate the impact of these lattice artifacts on further observables like magnetic polarizabilities and discuss the magnetic field-induced mixing between $\rho$-mesons and pions. We also derive Ward-Takashi identities for QCD+QED both in the continuum formulation and for (order $a$-improved) Wilson fermions.


I. INTRODUCTION
Background magnetic fields have a decisive impact on the physics of quarks and gluons and offer a wide range of applications. Strong magnetic fields appear in noncentral heavy-ion collisions [1,2], inside magnetars [3], and might have been generated during the evolution of the early Universe [4]. An important characteristic of the magnetized QCD medium -relevant for each of the above settings -is the effect of the magnetic field on the spectrum of the theory. It has been speculated [5,6] that long-lived magnetic fields might affect hadronization in heavy-ion collisions. In strongly magnetized neutron stars the magnetic properties of QCD matter, predominantly those of the light hadrons, may have a significant impact on the mass-radius relation [7].
The direct influence of the external field on the phase structure of QCD has received considerable attention recently; see for instance the review [8]. In particular, lattice simulations have revealed that the transition temperature of the chiral restoration/deconfinement crossover monotonously decreases as the magnetic field grows [9,10]; see also Refs. [11,12]. Thinking in terms of a hadronic model, in which the temperature needed to resolve the quark structure of a hadron is of the order of the hadron mass, it is plausible that the magnetic field dependence of the masses of the lightest hadrons do play a relevant role for the structure of the phase diagram. Specifically, pions are expected to be dominant in this respect. It is well known from chiral perturbation theory [13] that charged and neutral pions respond oppositely to the magnetic field B = B . While lattice results employing staggered [9] and quenched unimproved Wilson [14] quarks agree that the charged pion mass m π ± is increased by the presence of a magnetic field, there is a discrepancy in the literature about m π 0 (B) for intermediate magnetic fields. It monotonously decreases for quenched overlap quarks [15], while for quenched unimproved Wilson fermions, a turning point and subsequent increase was reported in Ref. [14]. 1 In light of the above discussion, the two results hint at opposite implications for the QCD phase diagram for strong magnetic fields.
Clarifying this question may also be relevant for yet another aspect of the phase structure of QCD with B > 0: It has been conjectured that strong magnetic fields could reduce the vector meson mass to zero and lead to the condensation of ρ-mesons [16]. They may also induce a transition to a superconducting phase [17]. 2 On the one hand, QCD inequalities can be used to show that massless ρ-mesons are only allowed if the (connected) neutral pion mass vanishes as well. The existence of a turning point in m π 0 (B), above which it grows with increasing B, would thus speak against this scenario. On the other hand, a study using quenched overlap quarks in two-color QCD, based on spatial vector correlators, has reported a critical magnetic field of around 0.9 GeV 2 , beyond which ρ-condensation was found to occur [20].
Resolving the above discrepancy regarding m π 0 (B) is thus clearly of high importance. In Ref. [21], we published preliminary results that explain the origin of the difference between the quenched overlap and Wilson results. We have found that Wilson quarks are subject to a B-dependent additive mass renormalization at finite lattice spacings. Although this effect disappears in the continuum limit, it has a drastic impact on the meson masses evaluated at nonzero lattice spacings. We also sketched a method to subtract these artifacts by determining a magnetic fielddependent line of constant physics for the Wilson bare mass parameter.
In this paper, we expand on the ideas introduced in Ref. [21] and discuss the improvement of scaling toward the continuum limit via the B-dependent additive mass renormalization in more detail. Performing this improvement reveals that the neutral pion mass monotonously decreases as B grows and, thus, resolves the existing disagreement between the quenched Wilson and overlap formulations. We also determine the ρ-meson mass and find that it remains nonzero in the whole range of magnetic fields studied, 0 < eB < 4 GeV 2 . For comparison, we also measure the connected neutral pion mass using dynamical staggered quarks on existing configurations from the study of Ref. [9]. 3 There are further interesting quantities that are related to the spectrum, including magnetic polarizabilities and magnetic moments. These describe the leading-order response of hadrons to the external field and, thus, involve derivatives at zero field; see, e.g., Refs. [15,22,23]. In contrast to observables at finite values of the field, these can also be defined for background electric fields, giving rise to electric polarizabilities. We will show that for Wilson quarks both types of polarizabilities suffer from enhanced lattice artifacts that can (and should) be eliminated using our improvement program.
We emphasize that our results are obtained in the quenched approximation, just as all studies in the literature so far concerning the light meson spectrum at finite B (except for the staggered results of Ref. [9], which, however, focused on pions only). The quenching of virtual quark loops induces a systematic effect that is hard to estimate. At B = 0 it is known that the results for the quenched and the dynamical spectrum agree within about 10% (see Refs. [24,25], for instance), indicating that the meson masses are dominantly determined by the properties of valence quarks. This is expected to remain true in the presence of external fields. For example the mass of the neutral pion, being a pseudo-Goldstone boson even for B > 0, depends primarily on the quark masses and not on gluonic properties. For further mesons one may expect a moderate quenching effect as long as the vacuum exhibits a structure similar to that at B = 0. Observables related to the bulk properties of the gluonic field (such as the static quark-antiquark potential considered in Ref. [26]), however, are expected to show stronger quenching effects, since the primary effect of the external field enters only indirectly via the sea quark loops which are not present in the quenched setup.
This paper is organized as follows. We start with a few remarks about the meson spectrum in continuum QCD in the presence of background magnetic fields in Sec. II. This is followed by Sec. III, where our lattice setup is described and our improvement scheme relevant for Wilson fermions is detailed. Section IV contains our results about the spectrum and the analysis of the dependence on quark masses and on the lattice spacing. In Sec. V we consider the impact of our improvement scheme for magnetic polarizabilities before we conclude in Sec. VI. Three Appendixes contain our results in the free case (Appendix A), the derivation of Ward-Takahashi identities for nonzero electromagnetic fields (Appendix B) and a perturbative check of the B-independence of the multiplicative mass renormalization constant in QCD (Appendix C).

II. MESON SPECTRUM IN THE CONTINUUM
A. Meson and quark masses for B > 0 We begin with a few general comments about the magnetic field dependence of the pion masses in the continuum. For small magnetic fields (and at zero temperature) we may approximate the charged pion as a pointlike free scalar particle with electric charge ±e. In that case, the quantum mechanically allowed energies are the Landau levels, where m π is the B = 0 pion mass and we have assumed B = B e z . Thus, the magnetic field-dependent mass (the lowest energy with n = p z = 0) reads The same argument for the neutral pion gives a magnetic field-independent mass, m π 0 (B) = m π . Thus, the neutral pion remains massless for nonzero magnetic fields if m π = 0, while the charged pions become massive for B ≠ 0. It turns out that the conclusion of increasing charged pion masses and an almost constant neutral pion mass remains valid if the interaction between pions is taken into account. This can be done consistently within chiral perturbation theory for nonzero magnetic fields [13,[27][28][29]. Indeed, the neutral pion remains a Goldstone boson even for B ≠ 0. In addition, away from the chiral limit, the Gell-Mann-Oakes-Renner relation still applies and is independent of B [27,28]. Here F π 0 is the decay constant of the neutral pion,ψψ is the (in our convention, positive) chiral condensate and m f denotes the (degenerate) up and down quark masses. F π0 andψψ above are defined in the chiral limit, m π0 → 0. Equation (3) is consistent with the notion that the quark mass does not depend on the magnetic field. The independence of the quark masses on the magnetic field also becomes obvious in a completely different limit: at very high temperatures. Here, quarks become quasifree due to asymptotic freedom. Then the quantum mechanical energy levels are the Landau levels for fermions with charge q, so that the lowest energy (n = p z = 0) is indeed magnetic fieldindependent. 4 The free case is of interest for yet another reason: In the strong magnetic field limit B → ∞, quarks and gluons decouple and pure fermionic observables approach their free-case values due to asymptotic freedom (see, for example, Refs. [10,32]). It is therefore instructive to compare our QCD results at large magnetic fields to the free case; see below.
In QCD, the quark mass is also subject to a multiplicative renormalization of the form m r f = Z m m f . Since Z m is related to the ultraviolet behavior of the theory, while the magnetic field is a physical infrared parameter, one expects that Z m is also independent of B. We demonstrate this by an explicit calculation using perturbation theory in Appendix C.
In the following, we will also investigate the effect of the magnetic field on ρ-mesons. Unlike for pions, for vector mesons, the magnetic field interacts directly with the spin of the particle. In the free case, the energies of a pointlike vector meson (with a gyromagnetic ratio g = 2 and charge q = ±e) are given by [19] where s z = 0, ±1 is the projection of the spin on the magnetic field axis. For s z = 1 one may expect the mass of the associated ρ-meson to vanish when eB = m 2 ρ . This has lead to the proposal that the ρ-mesons could condense [16] and initiated the investigation of the possibility of a superconducting vacuum for strong magnetic fields [17], as mentioned in the Introduction. However, QCD inequalities imply [14] that the correlation function in the vector channel is bounded from above by the (connected) neutral correlation function of the pion, so that Thus, a vanishing ρ-meson mass is only possible when m π u π d = 0. Note that the ρ-meson is a resonance in nature, rather than a stable particle, so the extraction of its resonant mass from lattice simulations demands a careful finite size analysis. In our study, however, the ρ-meson can be considered to be a stable particle, since we are working in the quenched approximation. When we consider the neutral pion at B = 0, where the SU V (2) symmetry is intact, the associated state is a mixture of neutral pions withūu anddd flavour contents, 4 Clearly, if QCD interactions are present, the levels will mix. Nevertheless, it turns out that the lowest Landau level can still be defined unambiguously and remains approximately B independent [30,31].
At B > 0, however, SU V (2) is broken explicitly by the different quark electric charges, q d = −q u 2 = −e 3, where e > 0 is the elementary charge. In this case, the neutral pion will be given by the more general relation where α(B) → 1 √ 2 and β(B) → 1 √ 2 for B → 0. Thus, determining quantities related to the neutral pion in principle involves the computation of the coefficients α and β for each value of the magnetic field. Note that for B > 0, the associated correlation matrix contains disconnected diagrams (these terms cancel at B = 0 due to isospin symmetry). In the present study, we neglect disconnected diagrams and the associated mixing and work with the individual π u and π d -states instead. This approximation is expected to become valid for strong magnetic fields, where the free case is approached (see the discussion above) and, thus, disconnected diagrams become negligible. Note that the so-defined connected pion states π u and π d are still useful as, for example, they enter in QCD inequalities [14]. We remark furthermore that neutral ρ-mesons are affected by the above issue in the same way.
The presence of a finite external magnetic field enables the mixing between the pion and the ρ-meson with s z = 0, Throughout the paper, superscripts of meson states denote electric charge, and subscripts denote the spin. For s z = ±1, the above mixing is forbidden due to the conservation of angular momentum. Consequently, pion states contribute to the spectral representation of correlation functions associated with s z = 0 ρ-mesons, and, in turn, these ρ-meson states also contribute to the spectral representation of pion correlation functions, so both correlation functions show the same leading-order exponential falloff. In fact, once this mixing is enabled, the physical mass eigenstates are mixtures of would-be π and ρ-states which can be written as where θ is the mixing angle and we denote the state with smaller energy by π ′ . In principle, the corresponding energies can be extracted from a correlation matrix using the operators for π ±,0 and ρ ±,0 0 , while θ will depend on the interpolators in use. On top of these two (and higher excited single-particle) states, there are also states containing more pions (and radial excitations) that may be involved in the mixing. However, their contribution is negligible if m ρ ′ < E 1 is satisfied, where E 1 is the energy of the lowest-lying multiparticle excitation.

A. Lattice setup
We consider two flavors of (quenched) unimproved Wilson fermions on the lattice with the fermion action Here a denotes the lattice spacing, ψ = (u, d) ⊺ is the Dirac field for flavor u and d, and M 0 is the matrix of bare quark masses, where 1 is the unit matrix in flavor space and τ 3 is the third Pauli matrix. γ µ denote the Euclidean Dirac matrices. r is the coefficient of the Wilson term which we will set to r = 1 in the simulations. The quarks couple to gluonic and electromagnetic fields via the link variables U µ (x), containing the gluonic links U G µ (x) ∈ SU(N ) and the electromagnetic link variables u µ (x), The latter are matrices in flavor space since the electromagnetic charges are different. We are mostly interested in QCD in the presence of an external magnetic field B = B e z . Such a magnetic field can be generated by a vector potential of the form where the parameter η cancels in physical observables and thus can be chosen at will. In the following results are obtained using the "symmetric gauge," η = 1 2, but we have explicitly checked that the results agree when we set η = 1. L µ = N µ a denotes the spatial extent in the µ-direction and N µ is the number of lattice points. The boundary δ-terms ensure periodic boundary conditions for the gauge potential and the constancy of the magnetic field in the x, y-plane [33]. Note, that in a finite volume, the magnetic field is quantized according to where the smallest quark charge, q d = −e 3 appears.

B. Quark mass renormalization
The term proportional to r in Eq. (11) is the well-known Wilson term, which is introduced to remove the doublers from the theory. This term breaks chiral symmetry explicitly so that the quark mass is subject to additive renormalization, besides the multiplicative renormalization that is already necessary in the continuum. From the arguments of Sec. II, we expect the multiplicative renormalization factor Z m to be independent of the external field -or, in other words, a B-independent renormalization scheme can be chosen (see Appendix C). Thus, at finite lattice spacing, the renormalized quark mass of flavor f is given by 5 where and d m is the coefficient of the leading-order B-dependent term. Note that, although d m corrects for a lattice artifact, it is not a Symanzik improvement coefficient. However, it can have a nontrivial value since our theory is not O(a 2 ) [and not even O(a)] improved. In Appendix A we give an intuitive explanation of why magnetic fields shift the quark mass, i.e., why m c;f becomes B dependent. The argument boils down to the fact that the Wilson term coincides with the action of a scalar particle -multiplied by the lattice spacing so that it vanishes for a → 0. In the presence of the magnetic field, the squared mass of this scalar particle is shifted by B according to Eq. (2). Eventually this induces the B-dependence in the additive mass renormalization. This can be checked explicitly in the free case (where d m = −1 2), for which we show results in Appendix A.
It is worth stressing that the presence of the additive quark mass renormalization is a lattice artifact and m c;f (a, B) → 0 for a → 0 for any magnetic field. For Wilson fermions it is customary to introduce the hopping parameter κ f and its critical value κ c;f , so that the B-dependent quark mass renormalization can be translated to the B-dependence of κ c,f , At finite lattice spacing the B-dependence of m c;f contaminates physical B-effects. Due to this additional lattice artifact, very fine lattices are necessary for a reliable continuum extrapolation. A better way is to get rid of the contamination by keeping the simulation on the line of constant physics (LCP), characterized by a constant renormalized quark mass m r f . This LCP depends on the magnetic field according to Eq. (16), and will be referred to as LCP(B) from now on. 6 Since Z m is, in our definition, B independent, it suffices to tune the hopping parameter κ f = κ f (a, B) such that the bare (multiplicatively unrenormalized) quark mass 5 Note that a similar equation holds for the case of the coupling to a dynamical electromagnetic field (e.g. Ref. [34]), m r f = Zm(a, e) m 0 f − m c;f (a, e) , where both renormalization factors depend on the bare electromagnetic coupling e. 6 We remark that, for unimproved Wilson fermions, lines of constant physics are, in general, only defined up to magnetic field-dependent or -independent O(a) lattice artifacts. Here, we choose to correct for the B-dependent artifacts originating from the additive mass renormalization, since these turn out to be numerically large.
remains constant when varying B. We will discuss the impact of this tuning on lattice artifacts in detail in Sec. IV E, where we compare the results obtained when keeping κ f constant while varying B to those obtained along the LCP(B). The B-dependence of the quark mass is also inherited by further quantities that are related to derivatives with respect to the magnetic field. Particular examples, namely magnetic moments and polarizabilities, will be discussed in Sec. V. Before we move on, let us collect some theoretical expectations concerning the behavior of m c;f (a, B). Parity symmetry ensures that m c;f (a, B) = m c;f (a, −B). Below we show numerical results in the quenched setup, where charged quarks are only included in the valence sector. In this case q u q d = 2 leads to C. Determination of the critical mass The key ingredient to fixing the LCP(B) is the B-dependence of κ c;f . To determine it, we can make use of the fact that the neutral pion mass vanishes in the chiral limit, even at finite B (cf. Sec. II). In our quenched setup this chiral extrapolation is complicated by the appearance of chiral logarithms in quenched chiral perturbation theory [35,36] (see also the right panel of Fig. 1 below). It is more convenient to use the current quark massm f obtained from the axial Ward-Takahashi identity (WI). In Appendix B we derive the axial and vector WIs in the continuum and for Wilson fermions for B > 0. While the charged WIs acquire additional B-dependent terms, 7 the neutral WIs are found to keep their B = 0 form, Here, (J A ) f is the (point-split) axial current operator with fermionic fields of flavors f (cf. Eq. (B24)), and P f the associated pseudoscalar density. The isosinglet WI suffers from the well-known axial anomaly, see Eq. (B17), which involves only gluonic contributions for background magnetic fields, and obtains contributions from disconnected diagrams, which we neglect. This corresponds to pretending that we have two different flavors carrying the same electric charge, in analogy to replacing π 0 by π u , thereby avoiding disconnected diagrams. In the following, we will thus determine κ c;f by requiring that the "neutral" current quark massm f , defined via Eq. (22), vanishes at κ c;f (B) for each value of B. The correlation functions have been projected to zero momentum in the z-direction. In addition, we have summed the correlation functions over the x, y-plane at the sink. For B = 0, this corresponds to projecting to zero momentum in the x-and y-directions. On the level of the electromagnetic vector potential, translational invariance is broken, and Fourier transformation is not well defined in the x, y-plane (see also Ref. [14]). We found that the summation enhances statistics but does not affect our results, as we have checked explicitly by comparing the results from summed correlation functions to those obtained without summation at a fixed position in the x, y-plane. For each value of a and B, we measuredm f for several values of κ and determined κ c;f (a, B) via a linear extrapolation of the form [cf.
Here C f is a proportionality constant, which relates the two estimates for the bare quark mass to each other and, thus, may contain (B-dependent) lattice artifacts. It can be written as where we have introduced the B = 0 proportionality factor (see, e.g., Ref. [40]) and d Z is the coefficient of the leading-order B-dependent lattice artifact in C f . 7 See also [37][38][39].     We note in passing that ideally one would like to use the charged axial and vector WIs for the tuning of sums and differences of quark masses. The corresponding identities are worked out including electromagnetic interactions [and also O(a)-improvement] in Appendix B. In this case no disconnected diagrams contribute. We have also looked at the charged WIs but could not detect plateaus before the signal got lost in the noise. A potential reason for this is that the correlation functions used to construct the charged WIs decay with a mass larger than the smallest mass in the spectrum (since m π + > m π 0 ) and thus suffer from the standard signal-to-noise problem. We plan to address this and investigate the charged WIs in more detail in the future.

D. Lattice setup and lines of constant physics
The simulations are performed in the quenched setup using the Wilson plaquette action [41] and heat bath and over-relaxation updates with a ratio of 1:4. The generation of configurations and the measurements of correlation functions have been carried out employing a modified version of CHROMA [42]. The run parameters are collected in Table I. We also list the value of the Sommer parameter r 0 a obtained from the interpolation provided in Ref. [43]. To convert to physical units, we use r 0 = 0.5 fm. The lattice sizes have been tuned to keep the physical volume fixed. This also allows us to achieve the same physical magnetic field -complying with the quantization condition (15) -for the different lattice spacings. The correlation functions have been computed using a Wuppertal smeared [44] source including (three dimensionally) APE smeared [45] link variables and a point sink. For the inversions, we have used the DFL-SAP-GCR solver introduced in Refs. [46,47].
To determine κ c,f , we have measuredm u andm d for several values of κ with pion masses between 1 GeV and about 400 MeV on each of these ensembles and for magnetic fields 0 < eB < 4 GeV 2 . For B = 0, the above setup satisfies m π f L ≳ 3 for all measurement points. However, as we will see below, the neutral pion mass decreases significantly with increasing magnetic field, so that finite size effects can potentially increase as B grows. To show that this is not the case, we study finite size effects in Sec. IV F.
We have extracted κ c;f (a, B) using a linear chiral extrapolation of the form given in Eq. (23). We show the results for am u (a, B; κ) together with the chiral extrapolations for several values of B obtained on the β = 6.0 ensemble in  The plot shows that a chiral extrapolation using the larger pion masses would lead to similar results. However, in the region of small pion masses we observe deviations from the straight line, possibly due to the presence of chiral logarithms. We show the behavior of C u (a, B) in Fig. 2. Toward the continuum limit, C u approaches unity as is clearly visible in the plot. In addition, for large values of B, we expect C u to approach the free fermion case for which again C u → 1. This is indicated by the increase in C u with increasing B. The increase seems to level off at around eB ≈ 2 GeV 2 , which could be a sign for growing lattice artifacts.
The results for κ c;u and κ c;d are tabulated in Table II. To estimate the systematic uncertainty of the linear ansatz within the present range of quark masses we have repeated the procedure by including a term quadratic in 1 κ−1 κ c;f in Eq. (23). The spread of the results from these two fits is given as the systematic uncertainty in Table II. To visualize the effect of B on m c;f , we introduce the parameter In this ratio and in our quenched setup, the leading B = 0 lattice artifacts and renormalization factors cancel, so this is the ideal tool to see to what extent the renormalized quark mass is affected by the shift in m c;f (a, B). Note, that if ∆m c;f (a, B) = −1, this means that the B-dependent part of the additive quark mass renormalization is negative and as large in magnitude as the quark mass associated with a pion with a mass of 415 MeV at B = 0.
The results for ∆m c;u and ∆m c;d are displayed for different values of B and the different ensembles in Fig. 3. In the left panel of the figure, we show ∆m c;u versus eB for the different lattice spacings [∆m c;d follows from these results and Eq. (21)]. The plot indicates that the mass shift is positive for low magnetic fields, in contrast to the negative shift in the free case [see Eq. (A1)]. However, the effects of higher-order terms in B soon become important and ∆m c;u (a, B) becomes negative and linear in B for strong magnetic fields -again in accordance with the expectation that we approach the free case in this limit. The figure also shows that the magnitude of ∆m c;u (a, B) decreases with decreasing lattice spacing -reflecting the fact that the magnetic field dependence of m c;f is a lattice artifact and thus vanishes in the continuum limit. The right panel displays the results in the (∆m u , ∆m d )-plane, and the individual data points indicate different values of B, starting from B = 0 (solid black square) and proceeding along the colored lines which are included to guide the eye.
The snail shell-alike shape of the curves mapped out in this parameter space demands some discussion. First, all curves start with a slope which is close to 2 (indicated by the dashed black line). This is expected due to the relation given in Eq. (21), which fixes the factor between the derivatives of m c;u and m c;d at B = 0. The direction of movement in the plane, however, is determined by the prefactor of the first term appearing in the expansion of m c;f (a, B), Eq. (17), with respect to B and we have seen above that it is positive. Deviations from the black line indicate the onset of higher-order terms in the expansion of m c;f (a, B). For large magnetic fields we expect to approach the free theory (see Sec. II), where the slope is again 2, cf. Appendix A. This behavior is also visible in Fig. 3. Next, let us discuss the change in shape for a → 0. While the slopes for B → 0 and B → ∞ are dictated by Eq. (21) and the free theory case, respectively, the extent of the curve depends on the lattice spacing. Remember that the B-dependence of m c;f (a, B) is a lattice artifact, so the coefficients of the expansion of m c;f (a, B) become smaller for a → 0, as can be seen in the left panel of Fig. 3. Consequently, the movement along the curve as a function of B becomes slower, and in the continuum limit, the curve will collapse to the starting point. To see this, compare the points of the curves marked with the black circles -these all correspond to the same physical magnetic field eB = 3.302 GeV 2 (N b = 10).
Using the results for κ c;u and κ c;d tabulated in Table II, we can now tune κ u and κ d along the LCP(B), so that the additively renormalized quark mass from Eq. (20) remains independent of B.  (4)(3) TABLE II. Results for κc;u and κ c;d obtained from the linear chiral extrapolation of the current quark masses. The first uncertainty is statistical, and the second is the estimate for the systematic uncertainty obtained from the spread to the result from an extrapolation with a polynomial of second order.

A. Computation of meson masses
The meson masses are extracted from the large Euclidean-time behavior of the correlation functions according to where O is either the pseudoscalar density (for the pion) or the vector current (for the ρ-meson), The charged mesons a = ± involve τ ± = (τ 1 ∓ iτ 2 ) 2, while our (connected) neutral mesons either have pureūu In the presence of the background magnetic field, ρ-mesons can be classified in terms of the projection of their spin onto the magnetic field, sB = s z B. In terms of the vector current operator (28), the spin eigenstates with s z = ±1 and 0 are given by (see, e.g., Ref. [48]) respectively. Thus, we can compute the relevant correlation functions for s z = ±1 via and for s z = 0 by We remark that charge conjugation symmetry and parity symmetry ensure that m π + = m π − and m ρ + 0 = m ρ − 0 for the modes with spin projection s z = 0 and m ρ + + = m ρ − − and m ρ + − = m ρ − + for the states with nonvanishing spin projection on the B-field axis.
To determine the fit range for the function (27), we have looked for the region of the minimal value of t included in the fit, t min , where χ 2 dof is close to unity and where the extracted meson mass does not depend on the particular choice of t min . For mesons with masses larger than the neutral pion mass m u m d (B) the correlation functions show the typical exponential decrease of the signal-to-noise ratio, so that the maximal possible value for t min decreases. In those cases we have less control over contaminations from excited states and we note that these effects should be further investigated in the future.

B. Meson masses from lines of constant physics
We have measured the meson masses along the LCP(B) using the tuned values of κ c;u and κ c;d described in Sec. III D. We work with quark masses corresponding to pion masses which are above 415 MeV at B = 0, so that m π L ≳ 3 at B = 0. In this case, finite size effects for the meson spectrum are expected to be small. A study of finite size effects at B ≠ 0 is presented in Sec. IV F. We start the discussion of the results with the quark mass dependence of the meson masses on the β = 6.0 ensemble, i.e., at the intermediate lattice spacing a = 0.093 fm.
In Fig. 4, we show the masses of neutral (left) and charged pions (right) versus B, normalized by the mass at B = 0. The neutral pion masses fall off strongly before they level off at intermediate values of B, decreasing monotonously throughout. At around eB ≈ 2 GeV 2 , the mass of the neutral pions has decreased to about 60 to 70% of its B = 0 value. The quark mass dependence of this behavior is rather mild but the effect is enhanced toward lighter quark masses. The charged pions show the increase with the magnetic field [cf. Eq. (2) for noninteracting pions], with the tendency to undershoot the free-case prediction for larger values of B. The same behavior has also been observed for dynamical staggered quarks [9] and perturbatively using one-loop pion-nucleon interactions [49].
Next, we investigate ρ-mesons. Of particular relevance is the component ρ + + (or equivalently the ρ − − ), which is the particle conjectured to condense [17]  described by the exponential form (27) even for the largest magnetic field eB m 2 ρ (0) ≈ 6.7. 8 The results for the masses of the ρ + + -mesons (left) and the ρ + − (right) are shown in Fig. 5 versus B, again normalized by the B = 0 mass. The plot indicates that the mass of the ρ + + decreases monotonously with B, but it starts to level off at around eB ≈ 2 GeV 2 and remains in the region of 50 to 70% of its B = 0 value. In terms of the inequality (6), this behavior is expected, given the flattening of m π u and m π d at large B. In fact, the ρ + + -meson mass does not saturate the bound from Eq. (6) but remains at a multiple of the pion mass larger than 1, as shown in Fig. 6 (left). The plot indicates that the ratio m ρ + + (B) m π 0 (B) remains constant within 20% with B and increases as the pion becomes lighter. This indicates that one moves further away from the saturation of the bound. The ρ + − masses show the increase expected from the free case, Eq. (5) with g = 2, but tend to undershoot the free-case curve for larger magnetic fields. Once more, this feature has only a very mild quark mass dependence.
In Fig. 6 (right), we show the masses of the neutral ρ u ± -meson with nonzero s z . Since the neutral meson has no magnetic moment, this mass remains unaffected by the magnetic field in the free case, and indeed we find that the mass is insensitive to the spin direction. However, we observe an almost linear rise with B, in this case almost completely independent of the quark mass. We suspect that disconnected diagrams might be important for this channel.
Next shown in Fig. 7. The plot indicates that in both cases only mild lattice artifacts are present. The result of the same analysis for the ρ + + and ρ + − -states is shown in Fig. 8 in the left and right panels, respectively. As for the pion masses, we observe only a mild lattice spacing dependence. Concerning the saturation of the bound (6), reducing the lattice spacing tends to drive the results away from saturation, as can be seen from Fig. 9 (left). This means that cutoff effects and quark mass effects act in the same manner in this respect. The mass of the neutral ρ u ± -meson with nonzero s z , shown in Fig. 9 (right), is also mainly independent of the lattice spacing with the tendency to show noticeable effects for larger values of B.

C. Pion-ρ-meson mixing
We will now discuss the mixing between pions and ρ-mesons with spin projection s z = 0 along the magnetic field axis, as described in Sec. II B. Since we neglect disconnected diagrams, we will focus on the mixing between the charged meson pair π + and ρ + 0 . Note that taking into account this mixing is necessary for the determination of the mass of the heavier of the pair (the ρ + 0 ), but it does not affect the lighter state (the π + ). Due to the mixing described in Eq. (10), for B > 0, the naive operators for the states π + and ρ + 0 overlap. The orthogonal basis can be determined via the diagonalization of the correlator matrix  and the true masses are the ones associated with the temporal exponential decay of the eigenvalues of this matrix in the t → ∞ limit. We remark that the off-diagonal elements need to be multiplied by ±i when we translate from the Euclidean to the Minkowskian Dirac matrices (for the spatial components these differ by a factor i so that γ E z = −iγ M z ). To extract the eigenvalues of the matrix, we use a generalized eigenvalue problem (see Ref. [50] and references therein) with starting temporal extents of t 0 = 3a, 4a and 5a for β = 5.845, 6.0 and 6.26, respectively. In the construction of the correlator matrix we have used smeared operators at the source and sink, so the correlator matrix ought to be symmetric. In practice, this only holds within the statistical uncertainty, and we found it beneficial to stabilize the computation of eigenvalues by replacing the off-diagonal elements by (the real part of) their average.
The smaller eigenvalue of the correlator matrix corresponds to the heavier state, i.e., to the ρ + 0 -meson. To extract the mass of the ρ + 0 -meson, we perform a fit similar to that in Eq. (27). We show the associated effective masses for the two eigenvalues for the β = 6.0 ensemble for eB = 0.660 and 2.641 GeV 2 in Fig. 10. Note, that the GEVP has been setup with t 0 = 4a and that the signal soon becomes lost in noise, so we have left out the results for the effective masses once uncertainties become overly large. For comparison, we have also plotted the result for π + from the analysis of the pion correlator without the use of the GEVP. The results indicate good agreement for the mass of the ground state in this channel. The final masses for different values of the quark mass are shown in Fig. 11 (left), normalized to the ρ-meson mass at B = 0. We see that the mass of the ρ + 0 increases and is in agreement with the free case prediction, the colored lines, even though it shows the tendency to overshoot the prediction for intermediate values of B and to undershoot for B ≳ 3 GeV 2 . This agreement lends support to our method and suggests that the contamination from other excited states is suppressed. The dependence of the mass on the lattice spacing is shown  in Fig. 11 (right), revealing that, as before, lattice artifacts are smaller than our statistical errors. Note once more, that we are working in the quenched approximation, where sea quarks are absent and the ρ-meson is a stable particle. Unlike in the case with dynamical fermions (e.g., Refs. [51,52]), we thus do not need to consider multipion states in our analysis.

D. Results in the continuum limit
We proceed by performing the continuum extrapolation of the meson masses. As in the previous section we will work at a fixed B = 0 pion mass of 415 MeV. Since we employ unimproved Wilson fermions, meson masses should show lattice artifacts of O(a). Consequently, our continuum extrapolation is carried out linearly in a and we extrapolate the meson masses for each value of B individually. The results for this continuum extrapolation are shown for two representative cases [the neutral (uu) pion masses and the ρ + − masses] in Fig. 12. The linear extrapolation works well for most of the cases, giving χ 2 dof values in the region between 0.5 and 2. However, we cannot exclude that our data still receive significant contributions by terms of O(a 2 ). To estimate the associated systematic uncertainty we have performed another linear continuum extrapolation including only the points at the two smallest lattice spacings and use the difference between the two extrapolations as the systematic uncertainty. In Fig. 12 the resulting uncertainty, including the systematic part, is represented by the error bars of the points on the a = 0 intersect. The continuum extrapolations at high magnetic field suggest that O(a 2 )-effects become more pronounced as B grows. In this region the systematic errors are potentially underestimated and should be checked in future simulations on finer lattices.
We show the results for the continuum extrapolated meson masses in Fig. 13. The results basically confirm what we already found at finite lattice spacing in the previous section. The masses of the neutral pions (withūu anddd flavor content) decrease down to about 60%-70% of their B = 0 value at B ≳ 2 GeV 2 , while the masses of the charged pions increase, in agreement with the energies in the free case, Eq. (1). The data for m π u can be qualitatively described by a curve of the form Performing a fit to the data, including also the data for m π d , with B → B 2 following Eq. (21), we obtain a 1 = 3.2(8) GeV −2 and a 2 = 4.8(1.2) GeV −2 . The resulting curve is shown in Fig. 13 (top), too. The mass of the ρ + + , similarly to that of the neutral pion, decreases down to about 60% of its B = 0 value, where it starts to level off. In fact, the ρ + + mass remains at twice the mass of the π u . The masses of the other ρ-mesons with s z = ±1 increase, even those of the ρ u d ± -mesons, the B-dependence of which is, as for the neutral pions, only an indirect effect of B, due to the nonvanishing polarizability of the meson.

E. Comparison to trajectories with B-independent κ
In the previous sections, we have discussed the results for meson masses obtained along the LCP(B), where the bare quark mass has been tuned in order for the renormalized quark mass to remain constant with varying B. In the literature, this tuning has not been considered so far, and we will now show the problems that can arise without it. Of course, results obtained on the LCP(0) should also give the correct continuum limit; however, the lattice artifacts in this case strongly depend on B, as we will demonstrate. In this section we will focus on the results for the neutral pion with uu flavor content, since these are the ones which are affected most by changes in the quark mass. However, other meson masses are influenced similarly by this effect.
In Fig. 14, we show the results for the masses obtained with constant κ for different values of B compared to those calculated along the LCP(B) for different lattice spacings, corresponding to a B = 0 pion mass of 415 MeV. The plot indicates that the results with constant κ suffer from enormous lattice artifacts. These act in a way in which the meson masses are initially below the continuum result, while they overshoot it for strong B. This tendency originates from the nonmonotonous behavior shown in Fig. 3, where the critical mass initially becomes bigger, before it starts to decrease with B. To investigate the different lattice artifacts in more detail, we plot the masses against the lattice spacing for some representative cases in Fig. 15. The plots show that the two sets of results seem to converge toward each other as a → 0. However, in most cases, a controlled (linear) continuum extrapolation is not possible for the results at constant κ. The exception is for eB = 0.660 GeV 2 , where both linear extrapolations point to the same continuum limit. Summarizing, the results for constant κ show lattice artifacts that are much larger than those of the runs along the LCP(B), which is particularly true for large values of B. In Fig. 16, we compare the ρ + + -meson mass results obtained from constant κ and from the LCP(B), revealing the same tendency as observed above for the pions.

F. Finite size effects
In the previous sections, we have seen that the neutral pion mass decreases down to 60% of its B = 0 value with increasing B. Consequently, finite size effects will potentially increase, since m π u (B)L becomes smaller. Indeed, for the strongest magnetic field, this combination drops to a value of m π u (B)L ≈ 2, indicating that finite size effects can potentially be sizable. To investigate this systematic uncertainty, we have generated configurations on a 48×24 3 lattice at β = 6.0, providing a second, larger volume for this β value compared to the 48 × 16 3 lattice listed in Table I, which has been used in the previous sections. On these configurations, we have computed the meson masses for eB = 1.320, 2.641, and 3.962 GeV 2 -these magnetic fields can be realized with integer flux quanta N b on both volumes; see Eq. (15). For the computations, we have used the κ values obtained from the LCP(B) determined on the smaller volume. The results for the two different volumes are shown in Fig. 17 for the neutral pion and the ρ + + , which are the two most important mesons for our analysis. The plot indicates that finite size effects are small. For both mesons, the results change by around 2%-3% which is well below the typical uncertainties of more than 10% of our continuum extrapolated results. The same is true for the other meson masses. Thus, we conclude that increasing B does not enhance finite size effects and that these at present can be neglected in comparison to other systematics.

A. Background field method and additive quark mass renormalization
Further interesting features of hadrons are their magnetic moments and the electric and magnetic polarizabilities. The latter describe the indirect response of a bound state to external electromagnetic fields. The background field method gives access to both the magnetic [53] and the electric (see Refs. [54,55] for instance) properties. In fact, in the past few years, a number of groups have started to measure these quantities in lattice QCD using Wilson fermions; see Refs. [23,[56][57][58][59][60][61][62][63][64][65][66]. So far, the change of the additive renormalization of the quark mass has been neglected in these studies, each of which were carried out at a single lattice spacing. Before investigating the magnitude of this effect, let us first discuss the implications that the unwanted change in the renormalized quark mass in the presence of an external electric or magnetic field may have in this respect.
The energy of a relativistic particle (hadron) with mass m and charge q in a magnetic field aligned in z-direction is to leading order given by (e.g., Ref. [67]) where g H is the g-factor of the hadron H and β H is the polarizability, which is absent in the free field cases discussed in Sec. II. Note that our convention for the polarizabilities is in agreement with the ones from [23,65,67], which use the non-relativistic limit, and it includes a contribution from the tensor polarizability for a particle with spin [23,65]. The ellipses stand for higher-order terms in eB . The magnetic moment µ of the particle is related to its g-factor in the nonrelativistic limit. The g-factor and the polarizability represent the leading-and next-to-leading-order responses of the bound state with respect to an external magnetic field. In particular, and β H = 1 8πm Here, the derivative ∂ is defined as the derivative with respect to eB while keeping all other renormalized parameters fixed. When we measure these with Wilson fermions via the background field method, however, we obtain the full derivatives, which receive contributions from the implicit B-dependence of the quark mass, so that Similarly, additional terms emerge for the second derivative as well. Analogous relations hold also for measurements of the electric polarizability. In Eq. (36), the sum is over the quark flavors f present in the measurement. In the quenched setup the sum includes only contributions from valence quarks. For dynamical quarks exposed to the external field, however, the sum also includes effects from sea quarks. The additional terms including derivatives of the quark masses with respect to the external field are lattice artifacts and thus vanish in the continuum limit. They are not present for actions where the quark mass is protected by a remnant of chiral symmetry. For Wilson fermions, however, those terms are present and contaminate the results for magnetic moments and polarizabilities. These contaminations can be removed by tuning the quark mass along LCP(B)s as discussed above, resulting in a cancellation of the derivatives of the quark masses with respect to the external field. Even when the tuning is only approximate, it leads to a strong reduction of the additional terms in Eq. (36) and, consequently, a reduction of lattice artifacts. We would also like to emphasize that this effect will always be present (and of the same size), no matter how small the external field is since polarizabilities are derivatives with respect to the field, Eq. (35).

B. Magnetic moments and polarizabilities of mesons and the impact of LCP(B)s
We proceed by investigating the strength of the effect mentioned above for polarizabilities and g-factors of mesons. Our setup is not optimal for the extraction of these quantities, since the volume in our simulations is rather small and, thus, the smallest available external field following Eq. (15) is relatively large. In fact, we could only use the lowest two values of the external field for the extraction of g H and β H . Obviously, this practice does not lead to a precision determination of these quantities, but it is sufficient for the purpose of demonstrating the improvement achieved.
Let us start with the polarizability of the neutral and the charged pion. Concerning the neutral pion, we can only determine the polarizability of pions withūu anddd flavor content individually. The results versus the lattice spacing are shown in Fig. 18 (left), both the ones obtained from LCP(B)s (filled symbols) and the ones obtained by keeping κ constant (open symbols). At finite lattice spacing, there is a visible difference between the two sets of results, and lattice artifacts appear to be larger for the results from constant κ, as expected. We did not attempt a continuum extrapolation, but it is clear that a naive linear continuum extrapolation leads to different results for the improved [LCP(B)] case and the unimproved one. However, the data indicate that a linear continuum extrapolation is not valid in the regions of lattice spacings at our disposal. A careful extrapolation to the continuum should, eventually, lead to the same result for both cases. Comparing our results from LCP(B)s to those of Ref. [67], we see that they are in the same ballpark, in particular, those for the pion with dd flavor content. Concerning the charged pion, the rather large uncertainties for the masses at finite B lead to large uncertainties for the polarizability, so we cannot draw any conclusions from them and we will not discuss these results here.
Moving on to the ρ-mesons, the polarizability for the ρ u d with s z = ±1 and the ρ ± 0 can be extracted analogously to the ones for the neutral and the charged pions, respectively. However, once more the uncertainties are too large to draw any definite conclusion, so we exclude these results from the discussion as well. The situation is somewhat better for the g-factor of ρ + with s z = ±1. In this case the energies in Eq. (34) allow for the separation of the effects from the g-factor and the polarizability by building the combinations  Following Eq. (34), ∆ is proportional to g, while Σ only contains the polarizability β ρ + as a free parameter up to O( eB 2 ). Note that in the quantity ∆-effects from the additive quark mass renormalization will cancel, since they appear with equivalent prefactors in E 2 Consequently, we expect to obtain similar results for g ρ + from LCP(B)s and with constant κ. The results are shown in Fig. 18 (right). Indeed, the two sets of results show good agreement within uncertainties. The results indicate that g ρ + is close to the free-case value g = 2, which is in qualitative agreement with previous lattice computations [68][69][70] and findings in chiral perturbation theory [71].
Since the relative size of the change in quark mass depends on the bare quark mass we have at B = 0 (the change in m c;f is independent of m f ), we expect the impact of the additional lattice artifacts to become increasingly important for smaller quark masses. We test this intuition by looking at the polarizabilities of the neutral pions for different values of m π (B = 0) on the lattice with β = 6.0. The results are shown in Fig. 19. The increasing impact when going to smaller quark masses is clearly visible in the data. We note that lattice artifacts of this type could potentially also be responsible for the pion mass dependence of the electrical polarizability of neutral pions obtained using the connected part of the correlation function only [59,63,66]. It has been found that the polarizability becomes negative for physical quark masses, which disagrees with expectations from chiral perturbation theory [72][73][74].

VI. CONCLUSIONS
In this paper, we investigated the meson spectrum in QCD at zero temperature in the presence of background magnetic fields B with Wilson fermions in the quenched approximation. The new methods introduced in this paper allow us for the first time to perform the continuum limit for the Wilson spectrum at a finite value of the external field. The novelty of our approach is the introduction of a magnetic field-dependent improvement term that ensures that the renormalized quark mass remains independent of B. This requirement defines a magnetic field-dependent line of constant physics LCP(B) for the bare quark mass parameter. This B-dependent tuning of the bare quark mass is absent in the continuum (as we discuss in Appendix C) and only concerns fermion discretizations that suffer from an additive renormalization, such as Wilson quarks, where this tuning is beneficial already in the free case (see Appendix A). Note that a similar tuning along LCP(B)s should also be carried out for nonuniform external fields; see the discussion in Appendix A. The LCP(B) was determined using the lattice Ward-Takahashi identities for nonvanishing electromagnetic fields, which we derived in Appendix B.
We emphasize that the improvement only differs from the naive approach by lattice artifacts, which, however, may be large. In particular, we demonstrated that without the improvement meson masses suffer from enormous discretization effects for strong magnetic fields. Besides the impact for the strong field limit, we also considered the effect of the improvement for derivatives with respect to the magnetic field at B = 0 -i.e., magnetic polarizabilities. Also, here we found that our approach suppresses lattice artifacts and enables a flatter, controlled continuum extrapolation.
Our most important results about the spectrum involve the mass of the neutral (connected) pion, which was found to decrease monotonously as the magnetic field grows, and the mass of the lightest charged ρ-meson, which remains nonzero for the whole range of magnetic fields that we considered. The latter might be relevant for the ρ-meson condensation predicted to set in for strong magnetic fields [16,19]. Nevertheless, we mention again that our quenched analysis of the ρ-meson mass in constant background magnetic fields cannot capture the subtle details of the superconducting vacuum [17,18]. Indeed, there are examples in the literature which show that the interplay between sea and valence quarks can lead to unexpected effects [75]. Eventually, the study of the meson spectrum should thus be repeated including sea quarks. Note, however, that this demands the generation of new configurations for each value of the quark mass and the external field in the (minimal) N f = 1 + 1 setup at large values of N t , rendering the study with dynamical fermions at least 2 orders of magnitude more expensive than the present quenched study. In addition to the investigation of lattice artifacts and quark mass effects, we have also checked for finite size effects and found these to be negligible compared to the other uncertainties (see Sec. IV F).
Finally, we elaborate on the connection between the B-dependence of the mass of the lightest hadron and that of the QCD transition temperature. The latter was determined using continuum extrapolated dynamical staggered quarks with physical masses in Refs. [9,10]. Remember that according to our results the quantity m π u (B) m π (0) has only mild quark mass dependence, cf. Fig. 4. It is thus sensible to compare our continuum extrapolated results to the staggered results using physical pion masses. Employing the same definition involving connected neutral pion correlators, the Wilson and staggered results 9 for m π u (B) m π (0) are shown in Fig. 20. For eB ≲ 1 GeV 2 , where both results are available, the comparison reveals small differences due to the different B = 0 pion masses and the inclusion of sea effects in the staggered case.
Regarding the QCD transition temperature, we considered the parameterization of the crossover transition temperature T c (B) defined using the inflection point of the average light quark condensate [10]. According to Fig. 20, the B-dependence of the curves looks qualitatively very similar -with an initial reduction followed by a saturation to around 60% to 70% of the B = 0 value. This supports the picture sketched in the Introduction, where we compared the finite temperature QCD transition to the "melting" of the lightest hadron state. Note that the quantitative difference for intermediate values of B could originate from effects due to the charged pions, which, in this region, are potentially light enough to have a direct influence on the transition temperature. We remark that, while the inclusion of charged sea quarks appears to have a marginal impact on the T = 0 pion masses, it is known to make a drastic difference for T c [8]. A possible explanation is that the transition temperature is driven by bulk effects encoded in the configurations (i.e., the action used to generate them), so an action leading to the correct pion mass is essential to observe the correct features of the transition. For meson masses, measurements are most affected by the properties of valence quarks, which reflects itself in the fact that the quenched spectrum reproduces the QCD spectrum up to about 10% accuracy (see, e.g., Refs. [24,25]). Nevertheless, it would be interesting to understand this subtle difference between measurements of properties of the phase transition and the hadron spectrum in more detail. other hand, for the charged correlation function, one of the quarks (that with the smaller absolute charge) is forced to have its magnetic moment antiparallel to B. The energy in this case should be E(B) = m u + m 2 d + 2 q d B . Our numerical results for the energies at fixed κ = 0.124 are shown in Fig. 22 (left). The energy of the neutral pion increases with the magnetic field, indicating the unphysical increase of the quark mass by the amount a 2 q f B 2. Tuning the hopping parameters along the trajectory (A1) instead, the neutral pion mass remains constant as it should; see Fig. 22 (right). Our results are also in good agreement with the expectation for the charged "pion." The shift in the bare quark mass basically follows from the magnetic field dependence of the lowest energy state of a noninteracting charged scalar particle. For a homogeneous background magnetic field, this dependence can be found analytically and is given by Eq. (1). We also checked numerically what happens if the magnetic background field is inhomogeneous. In particular, we considered an oscillatory field B(x) ∝ sin(2πkx L x ) with k ∈ Z and a half-half field B(x) ∝ (−1) x div (Lx 2) . In both cases the bare quark mass was observed to increase with growing magnetic fieldalthough not linearly, as for the homogeneous background, but quadratically in a 2 qB.

Appendix B: Axial Ward-Takahashi identities for QCD+QED
In this Appendix we will derive axial and vector WIs for QCD in the presence of electromagnetic interactions with two light quark flavors, i.e., for QCD+QED. We will start by (re)deriving the continuum WIs in Appendix B 1, which, to our knowledge, have so far only been discussed in Ref. [37], before we discuss the WIs for (unimproved) Wilson quarks on the lattice in Appendix B 2. Note that an account of the continuum WIs has already appeared in Ref. [21]. The WIs for O(a)-improved Wilson quarks can be derived along the same lines and will be discussed briefly in Appendix B 3. Note that the axial Ward identity for the domain wall fermion discretization has been derived in Ref. [38] and recently for twisted mass fermions at maximal twist in Ref. [39]. 10

Continuum Ward-Takahashi identities
In the continuum, the Euclidean fermion action including u and d quarks in the presence of QED interactions is given by with the covariant derivative M the matrix of bare quark masses in the continuum, given by Eq. (12); and the charge matrix We have denoted the gluon fields by G µ (x) and the electromagnetic (photon) field by A µ (x). Ward-Takahashi identities can be derived by varying the expectation value of some test operator O under transformations of the fermionic variables in the path integral of the form corresponding to infinitesimal transformations under SU V (2), and for infinitesimal transformations under SU A (2). Here, τ j is any of the three Pauli matrices for j = 1, 2, 3 and the unit matrix for j = 0, corresponding to isosinglet transformations. Using these transformations one obtains (with in the flavor nonsinglet case, where the second equality is valid as long as the operator O has no support at the point x. In the singlet case, one also has to consider the variation of the measure. Using Eqs. (B1) and (B2), we obtain The result for the first term is well known for the two types of transformations (see, e.g., Ref. [78]), where the curly brackets denote the anticommutator, so we only need to compute the additional terms stemming from the second term in Eq. (B7). Here, (J V ) j µ (x) is the (continuum) vector current from Eq. (28), including the Pauli matrix τ j 2, and (J A ) j µ (x) is the local axial vector current 11 (B10) 10 We thank Davide Giusti for drawing our attention to this. 11 We use the notation J V and J A instead of V and A in this section to avoid confusion with the electromagnetic vector potential Aµ and to distinguish between point-split and local currents in the lattice regularization.
Considering vector transformations, the additional terms are given by so together with we arrive at the continuum vector WI Here, ijk is the totally antisymmetric tensor for i, j, k = 1, 2, 3 and zero if either of i, j, k = 0, and we have introduced the scalar density With axial vector transformations, we obtain so we get for the continuum axial WI Here we have introduced the pseudoscalar densities H µν and F µν are the gluonic and electromagnetic field strength tensors, respectively; and Tr denotes the trace over flavor and color indices. The terms in the last line of Eq. (B17) are the ones associated with the Jacobian of the transformation from Eq. (B5) [79], known as the axial anomaly. For QCD+QED, there are two such terms, associated with the topological charge operators for the gluonic and the electromagnetic fields, respectively. Note that for external magnetic fields (and no electric fields), the electromagnetic anomaly is absent.

Ward-Takahashi identities for Wilson fermions
We will now derive similar identities for (unimproved) Wilson fermions using the fermion action from Eq. (11). In particular, we need to compute (∂ S F ) (∂α j X (x)) for the transformations from Eqs. (B8) and (B9). The basic idea to separate the resulting equations in parts which are already known and new parts is the same as in the continuum case. In practice, however, the separation is a bit less intuitive since the electromagnetic field enters the action via the electromagnetic link variables u µ (x); see Eq. (13). The separation can be done by bringing the Pauli matrices from the variation of theψ variables to the right-hand side of the link variables. 12 When the links u µ (x) are absent, that are a combination of WIs with τ 3 and 1, are unaffected by the presence of the QED interactions. 13 The terms denoted by (1) arise from the Wilson term (consequently, they are multiplied by a factor of r), and, as X j , they are proportional to a dimension-5 operator and vanish in the continuum. However, as X j , these terms also mix with the scalar (for the vector WI) and pseudoscalar (for the scalar WI) density, respectively, and thus lead to an additional additive quark mass renormalization. In particular, the presence of such a term in the vector WI indicates that the quark mass difference does not vanish when the bare quark masses coincide, meaning that u and d quark masses renormalizes with different additive terms. The terms denoted by (2) are operators of dimension 4 and survive the continuum limit. In fact, they exactly resemble the additional terms in the continuum WIs, Eqs. (B13) and (B17), proportional to A µ (x).
To compute current quark masses, terms of type (2) should be included explicitly in their definition, since these terms are also present in the continuum WI, while terms of type (1) should be left out (similar to X j ). For an example, consider the definition of the sum of u and d current quark masses via the axial WI evaluated for the matrix τ + from Sec. IV A. Using the pseudoscalar density P + from Eq. (B18) as the operator O in Eq. (B6), a suitable definition for the sum of current quark masses is given by where (2) refers to the insertion of the associated operators from Eq. (B22). Note that the electromagnetic link variables that need to be included in the point-split currents depend on the particular flavor matrix for which the WI is evaluated. As an example, let us once more consider the WIs for the matrix τ + . Writing we obtainψ

Ward-Takahashi identities for O(a)-improved Wilson fermions
For O(a)-improvement in QCD+QED, we need to include two Sheikholeslami-Wohlert [82] terms in the action, where H µν is (a suitable discretization of) the gluonic field strength tensor and C µν = QF µν is the electromagnetic one, multiplied by the electric charge matrix Q. The terms in the WIs which result from these additional terms depend on the particular choice of discretization. Most commonly used is the clover discretization of the field strength tensor, and similarly for C µν with U G replaced by u. Here, U G ±µ,±ν (x) [and similarly u ±µ,±ν (x) for C µν ] denotes the multiplication of links around a plaquette in the ±µ, ±ν-direction, starting from point x.

Appendix C: Multiplicative mass renormalization in magnetic fields
In this Appendix we demonstrate that the QCD mass renormalization constant is independent of the background magnetic field. On general grounds, it is expected that the ultraviolet divergent renormalization constants of the theory do not depend on physical parameters like the magnetic field. Here we show this using one-loop perturbation theory in continuum QCD. 14 The mass renormalization constant is obtained from the fermion self-energy, where the double line represents the free fermion propagator in the background magnetic field and the curly line represents the gluon propagator. The SU(3) generators are denoted by t a .
In the Schwinger proper time representation [84], the Feynman-gauge gluon propagator reads, where p 2 ⊥ = p 2 x + p 2 y and p 2 ∥ = p 2 z + p 2 t and throughout we will assume that qB > 0. In fact, Eq. (C3) is the Fourier transform of the translational invariant part of the propagator, and excludes the Schwinger phase. Inverse Fourier transforming Eq. (C1) and multiplying by the Schwinger phase gives the coordinate space representation of both sides [32,85]. Accordingly, the perpendicular p ⊥ are not real quantum numbers but referred to as pseudomomenta.
The propagator can be expanded in powers of the magnetic field. The consecutive orders read 14 For a similar calculation at nonzero temperature, see Ref. [83]. and can be inserted, order by order in the magnetic field, intoΣ. In color space, we only have the trivial factor ∑ a t a t a = C F 1 with C F = 4 3. The Dirac structure can be simplified with conventional techniques. The B = 0 contribution to the self-energy reads, after integration over p, This can be analytically integrated in t. The result at low s behaves like 1 s, such that an UV cutoff is necessary at the lower limit of integration 1 Λ 2 < s, and the integral diverges as log Λ. To obtain the mass renormalization constant, we need to evaluateΣ at the pole ofS 0 . This amounts to the replacement k → −im (compare to Ref. [79] and remember that we are working with Euclidean metric). We arrive at which is the well-known one-loop mass renormalization constant in QCD. Equation (C6) coincides with Schwinger's result for QED (replacing g → q and C F → 1). We proceed with the O(qB)-term. After integration over p, we get which can be integrated directly over t. The result behaves as s 0 at low s, and thus involves no ultraviolet divergences. Analytically integrating over s results in which is related to the anomalous magnetic moment of the quark [84]. Finally, we consider the O((qB) 2 )-term. Also, here, we are able to perform the p-integral, followed by the t-integral. The result is linear in s at low s and is thus again ultraviolet finite. In fact, since higher powers of the magnetic field always appear with higher powers of s [see Eq. (C4)], the ultraviolet region s ≈ 0 is getting more and more suppressed as we consider higher and higher orders in the magnetic field. Altogether, this shows that the B-dependent part of the quark self-energy is, to all orders in B, finite. Therefore, the mass renormalization constant contains no B-dependent divergences and one is free to use a B-independent renormalization scheme to define the renormalized mass.