Magnetic catalysis in the (2+1)-dimensional Gross-Neveu model

We study the Gross-Neveu model in $2+1$ dimensions in an external magnetic field $B$. We first summarize known mean-field results, obtained in the limit of large flavor number $N_\mathrm{f}$, before presenting lattice results using the overlap discretization to study one reducible fermion flavor, $N_\mathrm{f}=1$. Our findings indicate that the magnetic catalysis phenomenon, i.e., an increase of the chiral condensate with the magnetic field, persists beyond the mean-field limit for temperatures below the chiral phase transition and that the critical temperature grows with increasing magnetic field. This is in contrast to the situation in QCD, where the broken phase shrinks with increasing $B$ while the condensate exhibits a non-monotonic $B$-dependence close to the chiral crossover, and we comment on this discrepancy. We do not find any trace of inhomogeneous phases induced by the magnetic field.

We study the Gross-Neveu model in 2 + 1 dimensions in an external magnetic field B. We first summarize known mean-field results, obtained in the limit of large flavor number N f , before presenting lattice results using the overlap discretization to study one reducible fermion flavor, N f = 1. Our findings indicate that the magnetic catalysis phenomenon, i.e., an increase of the chiral condensate with the magnetic field, persists beyond the mean-field limit for temperatures below the chiral phase transition and that the critical temperature grows with increasing magnetic field. This is in contrast to the situation in QCD, where the broken phase shrinks with increasing B while the condensate exhibits a non-monotonic B-dependence close to the chiral crossover, and we comment on this discrepancy. We do not find any trace of inhomogeneous phases induced by the magnetic field.

I. INTRODUCTION
In recent years the study of strongly interacting quantum field theories exposed to external electromagnetic fields has received a significant amount of attention in the high-energy physics community. This is due to the fact that magnetic fields are believed to play an important role in a plethora of physical processes, such as heavy-ion collisions [1][2][3][4][5], the strong interactions within neutron stars [6][7][8][9], and at several stages of the early Universe [10][11][12][13] -see [14] for an extensive review.
Quantum Chromodynamics (QCD) is the theoretical framework underlying the strong interactions and -as such -describes the aforementioned phenomena. Since QCD cannot be studied using perturbation theory in the parameter regime of interest, one has to resort to non-perturbative methods, of which lattice quantum field theory is the most reliable one. Lattice simulations, however, suffer from the infamous complex-action problem, rendering the use of conventional Monte-Carlo methods impossible at finite density.
Needless to say, there are countless attempts aiming at circumventing the complex-action problem (see, e.g., [15]), but none of them has fully solved it within finitedensity QCD. In this work we employ a different approach altogether, using a low-dimensional toy model, the Gross-Neveu (GN) model [16], as an effective description of QCD. This is motivated by the fact that the GN model shares a number of important features with QCD, such as chiral symmetry and its spontaneous breakdown, or (in 3 dimensions or less) renormalizability [17,18].
It should be mentioned that there exist more realistic models, bearing a closer similarity to QCD than the one considered in this work, for instance models of the Nambu-Jona-Lasinio (NJL) [19] or quark-meson (see, e.g., [20]) type. Still, the simplicity of the GN model merits its * j.j.lenz@swansea.ac.uk † michael.mandl@uni-jena.de ‡ wipf@tpi.uni-jena.de use as a starting point for the search for a description of QCD using effective models, which may then be expanded upon. Furthermore, we mention that the GN model and variants thereof are also interesting from a condensed-matter perspective as they have been used successfully to describe certain one-dimensional and planar materials, such as polymers [21][22][23][24], graphene [25][26][27][28], and high-temperature superconductors [29][30][31]. One should, however, take some care in translating the results because the mapping of physical (non-relativistic) degrees of freedom to the fieldtheoretical description with emergent Lorentz invariance is not always straightforward -see [32] for one example.
Four-Fermi theories, including GN-type models, have been extensively studied in the literature with a variety of methods. A first -and often quite reasonable -approximation is given by mean-field treatments, which become exact at infinitely large flavor numbers N f and can be systematically improved by expanding in orders of 1/N f . Since we use the mean-field behavior as a guideline and as an important comparison for our lattice investigation, we summarize the most relevant known results for our system of interest in the following.
Early attempts to study the GN model in three spacetime dimensions including an external magnetic field were made in [33] and extended to finite temperature in [34,35]. It was found that the magnetic field is a strong catalyst of chiral symmetry breaking, enhancing the chiral condensate at both zero and non-zero temperature. This effect, termed magnetic catalysis, was explained in [36] to be caused by a dimensional reduction due to the applied field, similar to the effect of the Fermi surface in superconductivity -see also the reviews [37,38].
The goal of this work is to investigate whether the magnetic catalysis in the GN model is merely an artifact of the large -N f limit, where quantum fluctuations are suppressed, or is present in the full theory at finite flavor number as well. Work in this direction has already been performed using methods superior to the mean-field approximation, such as the functional renormalization group [39] or optimized perturbation theory [40], both supporting the presence of magnetic catalysis also at finite arXiv:2302.05279v2 [hep-lat] 13 Jun 2023 N f . However, we are not aware of any lattice studies of the GN model in an external electromagnetic field to be found in the literature at the time of writing, and we attempt to fill this gap. To this end, we perform extensive lattice simulations using overlap fermions at both zero and non-zero temperature and for vanishing chemical potential. The finite-density case will be studied in detail in an upcoming publication.
We note that -to the best of our knowledge -this is the first lattice Monte-Carlo simulation of GN-type models that uses overlap fermions (although theoretical considerations already exist in the literature [41]). Thus, we put considerable effort into working out the technical details and intricacies involved. However, we decided that they are better suited to be part of another planned publication with a more technical and analytical focus.
As we value the reproducibility of our results according to the FAIR Guiding Principles [42] (see [43] for a recent study about its status in our community), we provide access to our simulation results in [44]. Furthermore, the scripts used to perform our data analyses can be found in [45].
This work is structured as follows: In Sec. II we introduce the GN model in an external magnetic field and discuss its large -N f limit. Our lattice formalism is outlined in Sec. III, and our results are presented in Sec. IV. Finally, we discuss and critically analyze our findings and put them into perspective with respect to known QCD results in Sec. V.

II. ANALYTICAL RESULTS
The GN model in its most basic form is given by the Lagrangian [16] featuring N f flavors of fermionic fields, collected implicitly in the tuple ψ and self-interacting via a scalar-scalar channel with coupling constant g 2 . The sum over flavors is implied in (1).
In order to bring the model into a form amenable to our mean-field treatment as well as to lattice simulations, one introduces an auxiliary scalar field σ by means of a Hubbard-Stratonovich transformation. The semibosonized, but fully equivalent, theory then reads where we have furthermore coupled the fermions to an external vector field A µ and e denotes the elementary electric charge. For the remainder of this work we shall be concerned with a (2+1)-dimensional space-time and four-component spinor fields transforming in a reducible representation of the Dirac algebra [46]. This allows one to introduce a "fifth" 1 gamma matrix γ 5 , anti-commuting with all other γ µ . The Lagrangian (2) is then invariant under a discrete Z 2 chiral transformation: It is this chiral symmetry and its spontaneous breaking that will be our main concern in this work.
As an order parameter for chiral symmetry breaking we consider the fermion condensate ⟨ψψ⟩, which is proportional to the expectation value of the auxiliary field σ by means of a Dyson-Schwinger equation: In the limit of an infinite flavor number, the path integral defining the partition function of the model, is, after integrating out the fermions, reduced to the problem of minimizing the effective potential where we have assumed σ to be homogeneous in space and time, V denotes the space-time volume and D is the Dirac operator For A µ describing a constant and homogeneous (electro)magnetic field B and, without loss of generality, assuming B > 0, one finds the following effective potential density [48]: where σ 0 > 0 denotes the minimum of V eff at vanishing temperature and magnetic field, ζ H is the Hurwitz zeta function, β = 1/T denotes the inverse temperature, and the last term is a sum over Landau levels l with degeneracies d l = 2−δ l0 . Note that we work in the strong-coupling regime, where chiral symmetry is spontaneously broken at vanishing T and B, i.e., σ 0 ̸ = 0, which is not the case for weak couplings. Remarkably, the volume-dependence of V eff is contained only in the discretization of eB in a finite volume, see Eq. (19) below. For a derivation of Eq. (8), see App. A. The global minima ⟨σ⟩ of the effective potential for different temperatures and magnetic field strengths determine the mean-field phase structure of the GN model, which we show in Fig. 1. Evidently, chiral symmetry is spontaneously broken (i.e., ⟨σ⟩ ̸ = 0) for low temperatures and B = 0. The magnetic field then enhances this breaking even further, causing the chiral condensate to increase. This is the magnetic catalysis phenomenon mentioned in the Introduction. We furthermore observe that the critical temperature T c (B), beyond which chiral symmetry is restored (i.e., ⟨σ⟩ = 0), increases monotonically with B, and thus, the region of broken symmetry grows with the magnetic field.
We remark at this point that the magnetic-field-induced dimensional reduction down to one space-time dimension, found in [36,48,49] to be responsible for magnetic catalysis, is not in conflict with the no-go theorem prohibiting the existence of phases in one dimension [50] (not to be confused with the Coleman-Hohenberg-Mermin-Wagner theorem [51][52][53] preventing the spontaneous breaking of continuous symmetries in two dimensions). This is due to the fact that the chiral condensate itself is electrically neutral and, thus, unaffected by the dimensional reduction. For a similar argument in the U (2)-symmetric NJL model, see [48].
It is the main purpose of this work to shed light on the fate of the results presented in this section when going beyond the mean-field limit, i.e., when considering a finite number of fermionic flavors N f and lifting the restriction of homogeneity on σ.

A. Discretization
We intend to study the theory with Lagrangian (2) on a three-dimensional lattice Λ with N µ lattice points in the x µ -direction (µ = 0, 1, 2) and an isotropic lattice constant a. For the entirety of this work we shall always consider N 1 and N 2 to be equal, N 1 = N 2 =: N s , such that the physical lattice extent in each spatial direction is given by L = aN s . Furthermore, we introduce N t := N 0 to denote the number of lattice points in the (Euclidean) time direction, such that the inverse temperature reads β = aN t . We then denote the space-time volume as V = L 2 β. The bosonic field σ obeys periodic boundary conditions in all directions, while the fermions are periodic in space and anti-periodic in time.
The question of which lattice discretization to use for fermions is a non-trivial one. Studies of QCD in background magnetic fields mainly rely on the use of staggered fermions [54][55][56][57] (with a few authors employing overlap fermions as well [58]). However, it has become clear that staggered fermions can be problematic in asymptotically safe theories [59][60][61][62][63], of which (2) is an example. Moreover, since we are interested in studying chiral symmetry, we refrain from using Wilson fermions, and since we prefer to avoid the fermion doubling problem, we cannot use the naive discretization for N f < 8 either. Finally, even though in previous works [61,[64][65][66] the SLAC derivative [67,68] has proven to be the best-suited discretization for studying GN-like theories on the lattice, it fails when naively applied to theories with gauge symmetry [69]. As a matter of fact, it is not obvious how to properly formulate the GN model in a magnetic field with SLAC fermions in a gauge-invariant way in the first place. We nevertheless discuss this issue further and provide a more detailed comparison between different possible discretizations in App. B.
We are left with the choice of employing Ginsparg-Wilson fermions [70], which have ideal chiral properties but come with a significantly increased cost due to their non-ultralocality [71]. For our lattice studies we use Neuberger's formulation [72] of the overlap operator [73,74], Here, the kernel A is given by the Wilson operator D W with a negative mass parameter m = −1: 2 We remark that this expression does not make use of γ 5 and could thus be used in an irreducible representation of gamma matrices in (2 + 1) dimensions as well [75,76].
where the action of the covariant forward and backward difference operators on ψ(x) is defined as In (12),μ denotes the unit vector in the x µ -direction, and are U (1) link variables. Guided by the Lagrangian (2), where the Yukawa term σψψ would reduce to a fermionic mass term if σ was constant, one can introduce the scalar field into the overlap formalism by the definition [41] for the full Dirac operator. 3 For constant σ the second term in Eq. (14) is just a mass term in Ginsparg-Wilson language [78] (see [79] for a similar argument in the domain-wall formalism). By this definition one ensures that the identity (4), relating the expectation value of σ to the chiral condensate, is preserved 4 on the lattice, i.e., facilitating the numerical study of chiral symmetry breaking considerably. The full action of our lattice theory thus reads where summation over space-time and internal indices is implied. The discrete symmetry (3) of the continuum theory has an exact lattice counterpart in the overlap formalism, much like the case for theories with the more common U (1) chiral symmetry [80]. Namely, introducingγ 5 = γ 5 (1 − aD ov ), we find that the action (16) is invariant under by using the Ginsparg-Wilson relation [70] {D ov , γ 5 } = aD ov γ 5 D ov .
It should be noted that the additional symmetries in the continuum theory that arise due to ambiguity in the FIG. 2. Plaquette at position (x1, x2) in the spatial plane.
choice of the "fifth" gamma matrix (see footnote 1) can also be exactly translated to the lattice [81], but are not of interest in this work. From App. B we know that (massive) overlap fermions suffer from discretization effects that quantitatively change the chiral condensate in a theory of free fermions. Thus, one should investigate the interacting theory with a particular emphasis on its behavior towards the continuum limit to see if the discretization effects persist.

B. Magnetic field on the lattice
It is well known that the magnetic flux through a torus with an area L 2 , orthogonal to an applied magnetic field B, is necessarily quantized [82,83]. One finds the following quantization condition for the magnetic field: Let us now outline how to implement an external magnetic field perpendicular to the spatial plane using the gauge links (13) in our lattice formulation (16). In the continuum one could represent such a magnetic field by, e.g., the following choice of vector potential: On a lattice with periodic boundary conditions, however, this definition does not lead to a constant magnetic flux through every lattice plaquette P(x 1 , x 2 ) in the spatial plane at position (x 1 , x 2 ) -see Fig. 2 for the definition of such a plaquette and the integration path in Eq. (21). In fact, one finds i.e., the flux through the lattice boundary in the x 2direction is large and opposite to the flux through the bulk, such that the total magnetic flux through the lattice vanishes: The solution is to introduce correction terms in A µ on the lattice boundary in a way that shifts all the negative (assuming B > 0) flux to the single plaquette at the combined boundary x 1 = x 2 = L − a. This can be achieved by the following definition [84]: with A 0 set to zero. The flux through P(L − a, L − a) is now given by where we have used (19), and Φ P = a 2 B everywhere else. Since in our lattice formulation A µ only appears in exponentials due to (13), the only way Φ P contributes is via the plaquette terms as we have In this last expression the term proportional to 2π in (25) cancels out. We thus end up with a situation that is physically indistinguishable from one with a constant magnetic flux Φ P = a 2 B through every plaquette and a non-vanishing total flux as desired.
We therefore use the following definition of the U (1) gauge links U µ (x) in (13), entering the Wilson operator (11) via (12): We see that the compactness of the gauge links introduces a periodicity in the magnetic field and hence an effective upper bound for the flux quantum number b, i.e., In practice, one restricts b even further in order to avoid discretization artifacts [55,85] and we shall do the same in this work, performing simulations only up to b ≲ N 2 s /16.

C. Computational details
Our lattice setup of the GN model in a magnetic field, using the overlap Dirac operator (14), has a significant computational advantage compared to the use of overlap fermions in gauge theories. This is due to the fact that in our case the gauge links are not dynamical, depending only on the constant magnetic field. This allows for an exact computation of the massless overlap operator D ov in (9) that we perform once, at the beginning of a simulation. We then re-use D ov in every update step for the now straightforward computation of the full operator (14). Needless to say, computing the overlap operator exactly, i.e., without using approximations (see, e.g., [86]), would be unthinkable in realistic QCD simulations.
For this work we have performed simulations at various temperatures and magnetic fields using a standard rHMC algorithm. We change the temperature by varying N t at constant N s , and we study different lattice spacings by changing the coupling g 2 while simultaneously adjusting N s such that the physical lattice volume remains constant. We furthermore approach larger physical volumes by increasing N s at fixed g 2 . Finally, we mention that our theory does not suffer from a complex-action problem, as is shown in App. C.

D. Observables
As the order parameter for chiral symmetry breaking, the main observable of interest is the chiral condensate ⟨σ⟩ in (15). Assuming an ergodic simulation algorithm, however, this quantity will average to zero. This is because the effective potential of the GN model is known to exhibit two equivalent minima in the spontaneously broken phase, differing only in the sign of σ, hence leading to a cancellation between those minima. In order to avoid this cancellation, we thus use the quantity as an order parameter instead [87]. Here, the sum runs over the whole lattice and ⟨ · ⟩ denotes the Monte-Carlo average. While ⟨|σ|⟩ approaches ±⟨σ⟩ in the infinite-volume limit, one should keep in mind that on finite volumes ⟨|σ|⟩ will never be zero exactly, even when chiral symmetry is intact, which complicates the study of phase transitions.
For this reason, ⟨|σ|⟩ should -strictly speaking -not be referred to as an order parameter. However, for the sake of convenience we will still do so in the following.
In order to find the critical temperature T c , corresponding to the phase transition between the two respective regions of spontaneously broken and restored chiral sym-metry, we study the chiral susceptibility, defined as 5 as a function of T . Approaching a second-order phase transition, χ(T ) diverges rationally. This behavior is washed out by finite-volume corrections and we expect to find a sharp but smooth peak close to the transition temperature that monotonically grows, sharpens and moves towards the latter [88]. At this point we should mention that the introduction of an additional length scale and some form of imbalance in fermionic theories might induce spatial inhomogeneities in the system [89]. While this is most prominently observed in mean-field treatments at finite density [32,90], it could also apply to external magnetic fields. In fact, it is known that in 3+1 dimensions magnetic fields can favor inhomogeneous condensates at finite density when they would be disfavored at B = 0 [90][91][92][93]. This can be understood by recalling the dimensional reduction induced by the magnetic field [36] and the fact that inhomogeneous phases are more abundant in lower dimensions. Of course, our situation is qualitatively different because in our (2 + 1)-dimensional setup the dimensional reduction in a strong magnetic field leaves us with no spatial dimension at all (and we do not expect inhomogeneities in the temporal direction in equilibrium).
Since in 2+1 dimensions there is no conclusive evidence for the existence of inhomogeneous structures beyond mean-field (as compared to the (1 + 1)-dimensional case [64][65][66]94]) and there even exist some negative mean-field results [95][96][97], such inhomogeneities are not the focus of this work. Nonetheless, since the previous studies did not take into account the influence of magnetic fields, we also investigate whether an external magnetic field can induce inhomogeneities in 2 + 1 dimensions at zero density. To this end, we follow [64] by introducing the spatial correlation function As has been outlined in [64], this correlator should capture any inhomogeneities if they exist.

E. Scale setting
We set the scale via the order parameter at vanishing magnetic field and the lowest temperature considered, T 0 ≈ 0: 5 The factor of V is compensated by the use of space-time-averaged quantities in the expectation values such that χ is an intensive quantity, as it should be. We keep T 0 constant as we approach the infinite-volume (L 2 → ∞ at fixed a) and continuum (a → 0 at fixed L 2 ) limits, respectively. However, in order to ensure reasonably low scale-setting temperatures at an affordable computational cost, we consider two different T 0 corresponding to the two different limits.
For a detailed list of the parameters we have performed simulations for as well as their corresponding σ 0 and T 0 , we refer to Tab. I in App. D. In App. D we also give a brief description of how the error estimates presented in this work are obtained.

IV. RESULTS
In this section we report our results obtained in the GN model in 2 + 1 dimensions, using overlap fermions for one reducible fermionic flavor, N f = 1.

A. Consistency checks
As an important starting point, we test our discretization (14) and perform consistency checks with results in the existing literature. To this end, we show in Fig. 3 the dependence of the order parameter (31) on the coupling constant g 2 for increasing lattice volumes. The dashed blue line shows, exemplarily for the smallest lattice considered, the right-hand side of the Dyson-Schwinger equation (15) for comparison. This indicates that Eq. (15) is, indeed, fulfilled. The coupling strengths we use for the bulk of this work lie in the left half of Fig. 3.
In Fig. 3 we also show an extrapolation to the infinite volume, using the finite-size scaling law where α, γ and κ are constants, for the L-dependence of the order parameter for every value of the coupling. When a/g 2 takes values between 0.188 and 0.198 we find the offset α to be consistent with zero within errors in the infinite-volume limit, which indicates the presence of a phase transition. In this case, κ is related to the critical exponents β and ν of the order parameter and correlation length, respectively, via With this crude and naive method, we find β/ν = 0.93 ± 0.29 as a weighted average which -while not competitive in precision -is in quantitative agreement with pertinent results obtained with dedicated methods as, for example, collated in [98], β/ν = 0.62 . . . 0.85. Recovering this nonperturbative result is a strong indication that we are simulating the correct physics. From now on we shall always consider strong enough couplings such that chiral symmetry is spontaneously broken at T ≈ 0 = B, i.e., we work in the strong-coupling or super-critical regime as in Sec. II.

B. Vanishing magnetic field
Having established the correctness of our method, we now present results for the order parameter at vanishing magnetic field and non-zero temperature, which allows for a comparison with results in [87,99,100], at least on a qualitative level.
In Fig. 4 we show the T -dependence of ⟨|σ|⟩ for increasing physical volumes. We observe the expected spontaneous breaking of chiral symmetry at low temperatures, indicated by a non-vanishing order parameter, and a decrease of the condensate with increasing temperature, corresponding to the well-known picture of thermal fluctuations destroying long-range order and restoring chiral symmetry. Of course, as was mentioned above, ⟨|σ|⟩ cannot vanish exactly on finite volumes. What one can see, however, is that the phase transition becomes more pronounced as the volume increases, while the non-vanishing tail for high temperatures approaches lower and lower values.
In order to locate the phase transition we show in Fig. 5 the T -dependence of the chiral susceptibility (32) for different volumes. As expected, there is a pronounced peak at a critical temperature T = T c , which shifts slightly to lower temperatures as the volume is increased. For large enough volumes, where the peak becomes even more pronounced, we find T c /σ 0 ≈ 0.145.
We furthermore compute the Binder cumulant [101], as a function of T . The intersection of U L (T ) for different volumes provides us with another estimate for the critical temperature, T c /σ 0 ≈ 0.135. We take the interval between the two values as a rough estimate of the actual critical temperature. A direct comparison to the existing literature [87,99,100] is, unfortunately, not straightforward, as those works either employ higher flavor numbers or use different scale settings. The observations presented so far are consistent with the GN model approaching a second-order phase transition in T in the infinite-volume limit, as one would expect based on the large -N f analysis of Sec. II and as has been previously observed in [87,99,100].
Obviously, bosonic quantum fluctuations leave their mark on the system for flavor numbers as low as N f = 1, as can be seen by comparing the critical temperature quoted above with its large -N f value, T c /σ 0 = 1/2 ln(2) ≈ 0.72, the latter being significantly larger. This means that the broken phase shrinks when one departs from the meanfield limit by decreasing N f , which is not at all surprising given the tendency of quantum fluctuations to destroy any sort of long-range order. This phenomenon has also been observed in the earlier studies [87,99,100] and occurs in the (1 + 1)-dimensional model as well [64].
We remark that even the largest volume considered in this work is still comparatively small. Thus, one should not be tempted to draw quantitative conclusions about the precise location or the order of the chiral phase transition at N f = 1. The qualitative behavior, however, which is what we are ultimately interested in at B ̸ = 0, is as expected, which further builds up confidence in the chosen discretization.
C. Non-zero magnetic field

Temperatures close to zero
We now switch on an external magnetic field and first devote our attention to the lowest available temperatures. The B-dependence of the chiral condensate for various different lattice constants and volumes is shown in Figs. 6a and 6b, respectively. In all data the magnetic field is found to increase the chiral condensate which is in qualitative agreement with the large-N f expectation.
While the latter predicts quadratic growth for our scenario (and only linear growth in the sub-critical coupling regime), our data look rather linear but might still be compatible with a weak quadratic growth. This discrepancy could also come from discretization effects. Although Fig. 6a suggests that they may be small in the interacting theory, such a deviation would be the expected form of discretization artifacts in the non-interacting case as discussed in App. B. We find that such artifacts would systematically diminish the chiral condensate such that we are confident that our results are qualitatively correct even if discretization effects are larger than suggested by Fig. 6a.
Moreover, one observes curious non-monotonic behavior of ⟨|σ|⟩ with B, as the order parameter seems to assume a minimum at the lowest possible non-vanishing magnetic field, corresponding to b = 1 in Eq. (19), for all lattice spacings. For flux parameters larger than 1 the condensate then grows monotonically with B. This non-monotonicity, however, is a finite-size effect, as becomes clear by looking at the infinite-volume extrapolation shown in Fig. 6b, where b = 1 ceases to be a minimum of ⟨|σ|⟩ for the largest available volume (green curve). We note that the physical volume considered in Fig. 6a, which we keep approximately constant as we decrease the lattice spacing, corresponds to the smallest volume in Fig. 6b.
The largest magnetic fields we plot in Fig. 6 are determined by our requirement that b ≤ N 2 s /16. For larger magnetic fields we find unphysical saturation effects, the onset of which is already visible in the N s = 8 data of Fig. 6 (blue curves). We plan to present a more detailed discussion of these discretization artifacts and the aforementioned finite-size effects, as well as a thorough spectral analysis of the overlap operator for the GN model  in non-zero magnetic fields in a forthcoming publication.
We arrive at the conclusion that, on sufficiently large volumes and for temperatures close to zero, the magnetic field causes the order parameter to increase, thus enhancing the breaking of chiral symmetry, in accordance with the mean-field prediction of magnetic catalysis outlined in Sec. II. This is hardly surprising, given the effective one-dimensional dynamics induced by the magnetic field. In fact, as has been argued in [36], magnetic catalysis at zero temperature is a universal, i.e., model-independent feature in 2 + 1 dimensions, at least in the absence of gauge degrees of freedom [102].

Higher temperatures
Next, we study the combined influence of finite temperature and magnetic field on the order parameter. We show phase diagrams in the (B, T ) plane for various lattice sizes in Fig. 7.
Evidently, magnetic catalysis takes place not only for the lowest temperatures, but for all T below T c . We indicate the values of T c at B = 0 and N s = 16, determined above via χ and U L , respectively, by the gray bands. For higher temperatures the magnetic field ceases to have a noticeable effect on the order parameter. This is not unexpected, as in this region we only measure the (modulus of the) fluctuations of σ around zero due to our definition of ⟨|σ|⟩ in Eq. (31). The magnetic fields we restrict ourselves to (in order to avoid discretization effects) in our lattice simulations at fixed lattice spacing are quite small, eB/σ 2 0 ≲ 0.35. Hence, the results obtained in the large -N f approximation, shown in Fig. 1, suggest that one should not expect the broken region to grow in size all that much. Indeed, this expectation is confirmed by Fig. 7.
To investigate larger values of eB/σ 2 0 , we consider the (B, T ) phase diagram for the smallest available lattice spacing in Fig. 8. One observes that for strong enough magnetic fields the region of spontaneously broken chiral symmetry indeed starts to grow, as expected from Fig. 1. We roughly indicate this by the gray band, which shows the critical temperature T c , determined by the susceptibility (32), as a function of B. When T c cannot be determined unanimously we take the average of the two temperatures corresponding to the competing peaks instead and we do not show error bars for the resultingvery crude -estimate. Recall that finite-volume effects distort the behavior for weak magnetic fields. It would be interesting for future studies to consider even stronger magnetic fields in order to compare Figs. 1 and 8 on a more quantitative level. In conjunction with simulations at different flavor numbers, one could aim at finding a relation between the phase boundaries as N f is varied. In the simplest scenario, the critical temperature T c (B) could conceivably be related to its large -N f value by a mere N f -dependent scaling factor.

D. Search for inhomogeneous phases
Finally, we investigate the existence of inhomogeneous phases by studying the spatial correlator (33). Such a phase would likely occur at low temperatures and relatively strong magnetic fields, the former since thermal fluctuations will wash out any inhomogeneities and the latter since we know that the order parameter is homogeneous for vanishing magnetic field [95][96][97]. Fig. 9 shows the correlator C from Eq. (33) for a strong magnetic field along the two spatial coordinate axes and their diagonal. Each of them decays monotonically to a constant close to the contribution from the disconnected terms. In fact, we can showcase the rotational invariance here and no further structure is seen in other directions or for other parameters. We conclude that the assumption of spatial homogeneity is well justified in the accessible parameter range. Whether stronger magnetic fields could induce a spatially varying order parameter, especially in combination with a finite chemical potential, is a question for future studies.

V. DISCUSSION
We have investigated the (2 + 1)-dimensional Gross-Neveu model (2) exposed to an external magnetic field in the chiral limit using one reducible flavor of fermions and Neuberger's formulation (9) of the overlap operator. The auxiliary scalar field σ couples in a way so as to preserve the continuum Ward identity (4) relating the expectation value of σ to the chiral condensate.
Our results suggest that the magnetic catalysis phenomenon, i.e., an enhancement of the order parameter for chiral symmetry breaking with the magnetic field, persists for finite flavor numbers, in accordance with the phenomenological picture of the magnetic field reducing the number of spatial dimensions, thus promoting infrared dynamics. On small volumes, however, this effect is non-monotonic for weak magnetic fields. We also remark that our lattice formulation seems to suffer from strong discretization effects in a free-theory setup, while the interacting case appears less problematic.
We have furthermore investigated the fate of magnetic catalysis at finite temperature and found that it persists for all temperatures below the phase transition. Our findings are thus in qualitative agreement with mean-field [33][34][35] as well as beyond-mean-field [39,40] calculations. The phase of spontaneously broken chiral symmetry grows slightly for the strongest magnetic fields considered but shrinks overall in comparison to the large-N f limit.
It is important to stress that our results are very different from the well-known inverse magnetic catalysis effect, i.e., a decrease of the order parameter with B, that takes place in QCD at temperatures close to the chiral crossover [56,57]. In QCD, the critical temperature furthermore decreases with the magnetic field [103,104], which we do not observe either. We now comment on this issue.
In QCD the aforementioned effects are likely caused by a delicate interplay between quark and gluonic degrees of freedom [105], which our model, lacking the latter, cannot reproduce. Thus, one should not be tempted to interpret our results as new physics. Rather, we argue that the GN model is simply, and unsurprisingly, insufficient for a proper description of QCD once gluonic effects become important. While our results are thus in agreement with the expectation, we also stress that it was not entirely clear to us before this study to what extent the effect of quantum fluctuations might change the qualitative picture.
As was mentioned in the Introduction, we believe that our work serves as a starting point for the ultimate goal of studying QCD in background magnetic fields from the point of view of effective models beyond the meanfield limit and on the lattice. In the following we discuss ways to systematically improve the GN model in order to approach QCD. To this end, one should first consider models in 3 + 1 dimensions that have the same continuous chiral symmetry as QCD. One may then take gluonic interactions into account, for example, by coupling the fermions to the Polyakov loop [106]. Most importantly, the crucial back-reaction of magnetized quarks onto the gluonic distribution can be taken into account by introducing a suitable effective B-dependent coupling. This has been shown to reproduce the desired features of QCD in [107][108][109][110].
One could furthermore consider endowing the scalar fields in the NJL model with kinetic terms, thus enabling their interpretation as dynamical mesons and potentially add quartic mesonic self-interactions as well. The ensuing linear sigma model coupled to quarks (LSMq) has the added advantage that it is renormalizable in 3 + 1 dimensions, whereas the GN and NJL models are not. If one then incorporates the aforementioned magnetic-fielddependent couplings, while properly taking into account plasma screening effects, one also observes inverse magnetic catalysis as in QCD [111][112][113].
A proper understanding of such effective theories for QCD beyond the mean-field limit, e.g., from ab initio lattice simulations at finite numbers of quark flavors and colors is therefore certainly desirable. For reviews on the topic of reproducing features of QCD in magnetic fields using model theories and a more complete list of references, see [113][114][115].
We briefly comment on possible implications for condensed-matter systems that are described by four-Fermi theories. While in this work we are only concerned with the strong-coupling regime, in which chiral symmetry is broken at zero temperature and magnetic field, we believe that the qualitative predictions of mean-field studies should also remain valid for weak couplings. This would then imply that strong enough magnetic fields are indeed capable of generating a mass gap, providing fur-ther evidence [116,117] that magnetic catalysis could be responsible for the kink-like behavior of the thermal conductivity of superconducting cuprates exposed to a magnetic field observed in [118].
Finally, our results suggest that a small magnetic field does not seem to induce inhomogeneous phases in the GN model in 2 + 1 dimensions at zero density. A detailed study of the finite-density case is currently underway.
Our simulation results as well as the tools required to reproduce the figures shown in this work are available online [44,45]. This work would never have been possible without the great python ecosystem for scientific computing [119]. For our analyses, we explicitly imported the packages [120][121][122][123][124][125] but we are also grateful for creation and maintenance of all their dependencies.

OPEN ACCESS STATEMENT
For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any author accepted manuscript version arising.

DATA AVAILABILITY STATEMENT
Full data underlying this work are available in Ref. [44]. Fully automated analysis workflows are available in Ref. [45]. Raw data and the simulation code for generating the configurations are available upon request.
Appendix A: Derivation of the effective potential in the large -N f limit In this appendix we outline the calculation of the effective potential in Eq. (6), see also [126]. The main difficulty is, of course, the fermionic determinant, det(D). For the derivation we shall, in fact, consider the more general Dirac operator where we have also included a chemical potential µ and A µ is given in Eq. (20). In the following we assume B > 0 without loss of generality. For the computation of ln det D we use the zeta-function regularization method [127]: with the zeta-function of D 2 defined by where Γ(s) denotes the usual gamma function. The spectrum of D 2 is known and its eigenvalues read where ω n = π β (2n + 1) are the Matsubara frequencies (n ∈ Z), l ∈ N 0 denotes the Landau level index, and α ± 1 denotes the Zeeman splitting of energy levels of fermions with opposite spin due to the Pauli term in D 2 . The eigenvalues come with a degeneracy of 2 · V eB 2πβ , where the first factor of 2 comes from the use of a reducible representation of gamma matrices while the second factor is the standard Landau level degeneracy.
We are thus left with where we have already performed the sum over α and split up the summation over Landau levels into magnetic-fieldindependent terms (l = 0) and corrections due to B (l > 0). By performing a Poisson resummation in n and taking the integrals over t, a straightforward calculation leads to an expression for the zeta function, whose derivative with respect to s at s = 0 simplifies to where d l = 2 − δ l0 . After setting µ = 0 and inserting this expression into (A2) and (6), we obtain where we have replaced g 2 by the renormalized coupling g 2 R as dictated by the zeta function formalism. Finally, we introduce the minimum of the effective potential at vanishing temperature, density and magnetic field, σ 0 = −π/g 2 R , to recover (8).

Appendix B: Comparison of fermion discretizations
We discuss and compare three different fermion discretizations one could employ when attempting to study the GN model exposed to magnetic fields: naive, SLAC and overlap fermions. For the comparison we consider a theory of massive non-interacting fermions in an external magnetic field, characterized by the Lagrangian and compute the chiral condensate where the partition function Z is given by the fermion determinant, We have already computed ln det D in the continuum theory in App. A, allowing us to directly use the result (A6) by setting σ = m > 0. Thus, with Eq. (A2), the chiral condensate in the continuum at µ = 0 is given by the closed-form expression where ε l = m 2 + 2eBl .
We remark again that the volume-dependence only enters via the discretization of eB, Eq. (19). This means that if one were to naively take the limit eB → 0 in a continuous manner, one would simultaneously approach the infinitevolume limit.
To obtain the chiral condensate for vanishing magnetic field on a finite volume, one must repeat the calculation leading up to Eq. (B4), replacing the last term in Eq. (A4) by p 2 = p 2 1 + p 2 2 , with p i = 2π L n i and n i ∈ Z, and taking the sum over momenta in the place of Landau levels. Taking the four-fold degeneracy of the eigenvalues into account and repeating the steps outlined in App. A leads to the expression with for the chiral condensate on a finite volume and for B = 0. Let us now turn to the lattice computations. The basic ingredients for implementing external magnetic fields on the lattice are outlined in Sec. III B. For naive and overlap fermions we use the formalism developed there, mainly involving the U (1) gauge links in Eq. (29), which enter the naive Dirac operator, directly (∇ µ and ∇ * µ are defined in Eq. (12)) and the overlap operator via its Wilson kernel (10).
When using the SLAC derivative, however, one cannot use compact gauge variables in the form of group-valued lattice links connecting neighboring lattice sites because the derivative itself is non-local and thus involves all lattice points in a given direction. We therefore briefly discuss an alternative solution: In analogy to the continuum, we define the SLAC Dirac operator as where the SLAC derivative in position space is given by the Toeplitz matrix [128] ∂ SLAC if x µ ̸ = y µ and x ν = y ν for all ν ̸ = µ, and ∂ SLAC µ = 0 otherwise.
Obviously, the discretization (B9) is not gauge invariant. One could, however, attempt to treat the e / A term as a small perturbation if the magnetic field is not too large, such that (B9) still describes the correct physics approximately. 6 We do so in the following, reducing its numerical value as much as possible by employing the symmetric gauge with x 1,2 in the range −L 2 , L 2 . Problems will inevitably arise once the kinetic momentum p µ + eA µ crosses the boundary of the first Brillouin zone since there the SLAC derivative is discontinuous. This is the reason SLAC fermions are not used in gauge theories, and in our case such a crossing will occur for strong magnetic fields.
We note that the minimal coupling prescription used in (B9) makes the lattice boundary correction terms introduced in Eq. (24) obsolete, as they cannot be compensated for in the absence of compact periodic gauge variables. We have verified that their inclusion indeed gives worse results. Notice, however, that we are dealing with a different physical situation with SLAC fermions as now the total magnetic flux through the lattice vanishes, see Sec. III B.
Let us now compare the continuum result (B4) with the lattice chiral condensate, defined by where D latt stands for the operators D ov and D being defined in Eqs. (9) and (14), respectively. For naive fermions (B12) has to be divided by the number of doublers, i.e., 8 in 2+1 dimensions, in order to compare with continuum results. We show in Fig. 10 the change in the chiral condensate induced by the magnetic field, for the continuum result (where ∆⟨ψψ⟩ is obtained by subtracting (B6) from (B4)) and the three discretizations. One observes that the agreement with the continuum condensate is best for naive fermions. In an interacting theory, however, one cannot simulate the N f = 1 model with naive fermions by simply dividing by the number of doublers, which is the main reason we refrained from using the naive discretization in our study. The agreement for overlap fermions is very good for weak magnetic fields, in particular in the regime of eB/m 2 we investigate in our simulations. For stronger magnetic fields the qualitative behavior is still the same as for the continuum result, but the quantitative deviation (which appears to be quadratic in B) is substantial. We accredit this deviation to discretization artifacts, which for massive overlap fermions are worse (O(a)) than for naive fermions (O(a 2 )). One should therefore be cautious when interpreting our simulation results -while we do believe in their qualitative correctness, the absolute numbers could be systematically underestimated at large magnetic fields. In future studies one could employ an improvement program, such as the one suggested in [129], to reduce discretization effects.
For SLAC fermions, perhaps unsurprisingly, the agreement with the continuum result is rather poor, as the SLAC condensate does not even reproduce the qualitative features of the continuum one, for instance, the dip for the lowest allowed magnetic field. We mention a number of (ultimately futile) attempts we experimented with in order to improve the SLAC derivative in a magnetic field given in Eq. (B9).
First, we tried out different gauges instead of (B11), the latter leading to the best agreement, however. Next, we considered a physical situation where the magnetic field is constant and positive in one half of the lattice and constant and negative (with the same absolute value) in the other half. This avoids the need for introducing the lattice boundary terms in Eq. (24) entirely, which for the SLAC formulation were quite awkward in the first place. We then only considered the chiral condensate on a single lattice point x, lying in the center of the region with positive magnetic field. This was motivated by the intuition that at such a point the influence from the region with negative magnetic field should be negligible for large enough lattices. However, the agreement with continuum results we found was still poor. We conclude that more work is necessary if one aims at making SLAC fermions work for a background magnetic field.
where the σ µ can be chosen as the usual Pauli matrices. This decomposition makes clear how the reducible representation we use in this work is made up of the two inequivalent irreducible representations in three spacetime dimensions, σ µ and −σ µ . It is then straightforward to convince oneself that the overlap operator (14) also assumes a block form: where (i = 1, 2) and the irreducible components of the Wilson operator read (see Eq. (12) for the definitions of ∇ µ and ∇ * µ ) We emphasize that the diagonal elements D 1,2 in (C2) are precisely the expressions one would obtain for the overlap operator when working in one of the two irreducible representations. Hence, D decomposes in complete analogy to the continuum Dirac operator. Now, obviously, Furthermore, we note that the symmetric difference operator ∇ * µ + ∇ µ in (C6) is anti-Hermitian, while the discretized Laplacian ∇ * µ ∇ µ is Hermitian, such that D W,1 and D W,2 are Hermitian conjugates of one another. By using the spectral representation of the inverse square root in the definition of D ov,2 , one can then show that the same holds for D 1 and D 2 , such that, using Eq. (C7), i.e., there is no complex-action problem since the determinant is real and non-negative. We emphasize that the crucial ingredient for this proof was the use of a reducible representation of gamma matrices.

Appendix D: Parameters
In our simulations we generated O(10 3 ) − O(10 4 ) configurations per parameter set. We performed binned jackknife resamplings for our error analyses, making sure that each bin contained at least τ int configurations (but most commonly multiples thereof), where τ int refers to the integrated auto-correlation time corresponding to the order parameter ⟨|σ|⟩. We found τ int ≲ 50 in all cases.
We list the relevant parameters for which we have obtained simulation data, as well as the respective scales σ 0 and scale-setting temperatures T 0 , in Tab. I. Notice that the scale-setting temperatures are different between the infinite-volume and continuum limits, see Sec. III E. Since the errors in σ 0 are negligible, we do not quote them here and refrain from taking their influence on error propagation into account in the entirety of this work. I. Parameter sets used in the simulations. Here, Ns denotes the spatial lattice extent, assumed equal in both directions, Nt is the temporal extent, g 2 denotes the coupling constant in Eq. (16), b is the magnetic flux quantum number in Eq. (19) and T0 denotes the temperature at which we set the scale aσ0 in Eq. (34), which we quote in lattice units here. For the scale-setting data, the dots indicate steps of 0.005. As was explained in Sec. III E, we use different T0 for the infinite-volume and continuum extrapolations.