Neutron magnetic polarisability with Landau mode operators

The application of a uniform background magnetic field makes standard quark operators utilising gauge-covariant Gaussian smearing inefficient at isolating the ground state nucleon at nontrivial field strengths. In the absence of QCD interactions, Landau modes govern the quark energy levels. There is evidence that residual Landau mode effects remain when the strong interaction is turned on. Here we introduce novel quark operators constructed from the two-dimensional $U(1)$ Laplacian eigenmodes that describe the Landau levels of a charged particle on a periodic finite lattice. These eigenmode-projected quark operators provide enhanced precision for calculating nucleon energy shifts in a magnetic field. Using asymmetric source and sink operators, we are able to encapsulate the predominant effects of both the QCD and QED interactions in the interpolating fields for the neutron. The neutron magnetic polarizability is calculated using these techniques on the $32^3 \times 64$ dynamical QCD lattices provided by the PACS-CS Collaboration. In conjunction with a chiral effective-field theory analysis, we obtain a neutron magnetic polarizability of $\beta^n = 2.05(25)(19) \times 10^{-4}$ fm$^3$, where the numbers in parentheses describe statistical and systematic uncertainties.


I. INTRODUCTION
The study of the magnetic polarisability of the neutron is an area of ongoing experimental and theoretical interest. Measurement of this quantity remains challenging with considerable uncertainties [1][2][3] although improvement has been seen in recent years [4]. There is scope for lattice QCD to make important predictions in this area.
The approach used here to calculate this quantity on the lattice is the uniform background-field method [5][6][7]. A U (1) phase factor on the gauge links induces an external magnetic field across the entirety of the lattice. The external field causes an energy shift from which the magnetic polarisability can be determined by use of the energy-field relation [6][7][8][9][10][11] where m is the mass, and µ and β are the magnetic moment and magnetic polarisability respectively. Note that the |qe B|/2m term corresponds to the lowest Landau energy, and is only present for charged hadrons. There is in principle a tower of Landau levels, (2 n + 1)|qe B|/2m for n = 0, 1, 2, . . . [12]. At first glance, the method is simple; we can fit the linear and quadratic coefficients of the resulting energies as a function of field strength to extract the magnetic moment and polarisability [6,9]. However, baryon correlation functions suffer from a rapidly decaying signalto-noise problem [13]. This makes the extraction of the * ryan.bignell@adelaide.edu.au magnetic polarisability using standard nucleon interpolating fields challenging as it appears at second order in the energy expansion, as demonstrated by previous studies [7,[9][10][11].
The application of three-dimensional gauge covariant Gaussian smearing on the quark fields at the source and/or sink is highly effective at isolating the nucleon ground state in pure QCD calculations. However, the presence of a uniform magnetic field alters the physics, breaking three dimensional spatial symmetry, and introducing electromagnetic perturbations into the dynamics of the charged quarks.
Under a uniform magnetic field, in the absence of QCD interactions, each quark will have a Landau energy proportional to its charge. When QCD interactions are enabled the quarks will hadronise, such that (in the confining phase) the Landau energy corresponds to that of the composite particle. In particular, as the neutron has zero charge, the ddu quarks must combine such that the overall Landau energy vanishes.
It is clear that the QCD and magnetic interactions compete with each other in the confining phase. Indeed, there is evidence that residual Landau mode effects remain when the strong interaction is turned on [10,14]. Here we explore the idea of using quark operators on the lattice that capture both of these forces, choosing asymmetric source and sink operators to provide better overlap with the energy eigenstates of the neutron in a background magnetic field.
At the source, we capture the QCD dynamics by using spatial Gaussian-smearing, tuned to maximise overlap with the nucleon ground state at zero magnetic field strength. At the sink, we seek to encode the physics asso-ciated with the magnetic field by using a projection operator derived from the eigenmodes of the two-dimensional U (1) lattice Laplacian [15]. Calculations are performed at multiple quark masses in order to enable a chiral extrapolation to the physical regime.

II. BACKGROUND FIELD METHOD
To simulate a constant magnetic field along a single axis the background field method is used [5]. To derive this technique on the lattice, consider the continuum case. In the continuuum a minimal electromagnetic coupling is added to form the covariant derivative Here A µ is the electromagnetic four potential and qe the charge on the fermion field. The equivalent modification on the lattice is to multiply the usual gauge links by an exponential phase factor A uniform magnetic field along theẑ axis is obtained (in the continuum) using which does not uniquely specify the electromagnetic potential. The choice made over the interior of the lattice is A x = −B y. This gives a constant magnetic field of magnitude B in the +ẑ direction. In order to maintain the constant magnetic field across the edges of the lattice where periodic boundary conditions are in effect, we set A y = +BN y x along the boundary in theŷ dimension. This then induces a quantisation condition for the uniform magnetic field strength [10] qe B a 2 = 2 π k N x N y .
Here a is the lattice spacing, N x , N y the spatial dimensions of the lattice and k an integer specifying the field quanta in terms of the minimum field strength.
In this work the field quanta k is set in units of the charge of the down quark, i.e. q = −1/3. Hence a field with k d = 1 will be in the −ẑ direction and aligned with spin-down components.

III. QUARK OPERATORS
We explore the use of asymmetric source and sink operators in order to construct zero-momentum projected correlation functions which have greater overlap with the energy eigenstates of the neutron in a background magnetic field. This allows us to emulate the dominant QCD and magnetic effects separately.
We consider fixed boundary conditions in the time direction and place the source at N t /4 = 16.

A. Gaussian smeared source
A smeared source is used to provide a representation of the QCD interactions while the sink is used to capture the physics associated with the magnetic field. Several levels of source smearing are investigated at B = 0 in order to isolate the QCD nucleon ground state. For m π = 413 MeV 300, sweeps of standard Gaussiansmearing is optimal as illustrated in Figure 1. An identical process is followed at each of the quark masses producing optimal smearings of N sm = 150, 175, 300, 350 for masses m π = 702, 570, 413, 296 MeV respectively. B. Landau mode quark sink A charged scalar particle which sits in a uniform magnetic field will have an associated Landau energy which is proportional to its charge. In the non-relativistic approximation, the energy spectrum of a charged particle in a constant magnetic field along theẑ direction is equivalent to that of a harmonic oscillator, E n = (n + 1 2 )ω, where ω = |qeB|/m is the classical cyclotron frequency. In the continuum, each energy level is infinitely degenerate. The relativistic generalisation of the Landau energy levels for a fermion can be written as [12] where n describes the quantised energy level, α = ±1 is the spin polarisation eigenvalue, and p z is the momentum component in theẑ direction.
In an endeavour to capture the physics associated with the magnetic field, a quark level U (1) Landau projected sink is used. This sink aims to capture the effects of the magnetic field at the quark level. In the absence of QCD interactions, the charged quarks will each have an associated Landau level energy. The lattice Landau levels for a charged scalar particle correspond to the eigenmodes of the 2D Laplacian On the lattice, the degeneracy of the lattice Landau modes is finite and is dependent on the magnetic-field strength. This is in contrast to the infinite degeneracy of the continuum. In particular, the lowest Landau level on the lattice has a degeneracy equal to the magnetic flux quanta |k|.
The lowest Landau mode in the continuum takes a Gaussian form, ψ B (x, y) ∼ e −|qeB| (x 2 +y 2 )/4 . It has been noted elsewhere [9,16] that in a finite volume the periodicity of the lattice causes the wave function's form to be altered. We can calculate the eigenmodes of the 2D Laplacian in Eq. (7) and project at the quark level. Define a projection operator onto the lowest n eigenmodes | ψ i, B of the 2D Laplacian as A coordinate-space representation of this 2-dimensional projection operator is applied at the sink to the quark propagator, where n = |3 q f k d | modes for the lowest Landau level. The U (1) Laplacian is not QCD gauge covariant, and hence we fix the gluon field to Landau gauge and apply the appropriate gauge rotation to the quark propagator before projecting. However, as the hadronic correlation function (and ground state energy) is gauge invariant, using a gauge fixed sink operator can only effect the overlap with the ground state, which has the potential to improve the final precision of our result.

C. One dimensional spatial modulation
The eigenmodes of the two dimensional U (1) Laplacian have no dependence on the z coordinate. Using this freedom, we can apply a functional form to vary the spatial extent of the U (1) Landau projection in theẑ direction, an idea analagous to standard Gaussian smearing. We modulate the z dependence of the projected quark propagator with a normalised Gaussian, where the width parameter σ ≡ σ z controls the spatial extent in the z direction. After the U (1) Landau mode  projection has been applied at the sink to quark propagator as in Eq. (9), the gauge fixed propagator is then averaged over the z dimension using the modulation function as a weighting, We define the special case σ z = 0 to indicate that no z modulation is applied, which is equivalent to choosing φ σ=0 (z) = δ(z − z), such that S n,0 ≡ S n . Different spatial extents change the coupling to each of the energy eigenstates. The lowest lying level is dominant in the long Euclidean time limit. To determine which spatial extent provides the greatest overlap with the lowest lying energy level, many choices of σ z are investigated simultaneously [17].
The magnetic field orientation and neutron spin polarisation can be chosen independently to be in the positive or negative z direction. In order to efficiently extract the magnetic polarisability, combinations of correlation functions with differing magnetic field orientation and spin polarisation alignments are used to create spin and magnetic field aligned and anti-aligned correlation functions. These are the energies which will be examined in order to optimise the quark sink. The quark sink used is the one which has the longest plateau when fitting backwards in Euclidean time from where all of the correlators agree. This sink projected correlator that has converged the earliest is considered optimal. Figure 2 shows an example of this process for the m π = 411 MeV neutron and the largest magnetic field considered with |k| = 3. It is quite clear that all the sink projections agree by t = 29 and that σ z = 0, 1 both produce excellent plateaus. Of these two, σ z = 1 is preferred as it has a χ 2 dof which is closer to one. Note that the χ 2 dof is determined via a consideration of the full matrix of covariances between different time slices under consideration. Figure 3 shows the aligned energies for κ = 0.13770 in the smallest field strength. In this case there is no clear longest plateau and instead the χ 2 dof is the primary criteria. Here σ z = 1 is also preferred.
This process is undertaken for each combination of field strength and aligned or anti-aligned energies. The common result is that sink projections with σ z = 0, 1 are preferred. Hence these are the quark sinks considered when calculating polarisability energy shifts. These sink projections provide a good representation of the neutron ground state in a background magnetic field as can be seen by the rapid plateau behaviour in the energy of the neutron in Figure 4.
This result represents a significant advance in the determination of magnetic polarisabilities. For the first time clear plateaus are identified, a direct result of our consideration of Landau modes at the quark level.

A. Formalism
Recalling the energy-field relation of Eq. (1), we note that a combination of energies at different spin orientations and field strengths can be used to isolate the neutron magnetic polarisability β, noting that as q = 0 for the neutron the Landau energy term vanishes. Here the arrows denote the neutron spin polarisation along theẑ axis.  This method of isolating the polarisability term is valid, but in practice due to the cancellation of correlated fluctations on a common ensemble of lattice configurations it is much more effective to take ratios of appropriate spin-up (+s) and spin-down (−s) correlators. We can also average over both positive (+B) and negative (−B) magnetic field orientations to provide an improved unbiased estimator. Thus, we define the spin-field aligned correlator by and the spin-field anti-aligned correlator by The spin-field aligned and anti-aligned correlators, combined with the spin-averaged zero field correlator are used to form the ratio which upon taking the effective energy gives the desired energy shift Note that we define the magnetic field ±B to be that experienced by the nucleon, and is hence related to the down quark magnetic field by a factor of −3.
Any correlated fluctuations between the finite field strength and zero-field effective energies are significantly reduced prior to fitting by the zero-field denominator. The product of the spin-field aligned and anti-aligned correlators removes the energy contribution from the magnetic moment term. The use of aligned and anti-aligned correlators improves the statistics of the calculation by including the contributions from field and spin opposite pairs. This improved method is particularly important as the polarisability is at second order in B, and as such at these small field strengths its contribution to Eq. (1) is small. It is essential to have a precise determination of the polarisability energy shift. The efficiency of the Landau mode sink projection can be seen in Figure 5 where the energy shift for a standard, point sink is compared to a U (1) Landau mode sink projection; the latter is seen to display better plateau behaviour.

B. Simulation Details
In this work 2 + 1 flavour dynamical gauge configurations provided by the PACS-CS [18] group through the ILDG [19] are used. These have a clover fermion action and Iwasaki gauge action with a physical lattice spacing of a = 0.0907 (13). Four values of the light quark hopping parameter k ud = 0.13700, 0.13727, 0.13754, 0.13770 are considered, corresponding to pion masses of m π = 702, 570, 411, 296 MeV respectively. The lattice spacing for each mass was set using the Sommer scale with r 0 = 0.49 fm. The lattice volume is L 3 × T = 32 3 × 64 and the ensemble sizes are 399, 400, 449, 400 configurations respectively. Source locations were systematically varied in order to produce large distances between adjacent source locations. Starting from an initial source location at 0, 16 , shifts of 0, 16 were applied three times for a total of four source locations. A further set of four shifts starting at 16,8 , where only the time component increases by 16, were also applied. As such a total of eight sources were used for each configuration.  Correlation functions at four distinct magnetic field strengths are calculated. To do this propagators at ten non-zero field strengths, e B = ±0.087, ±0.174, ±0.261, ±0.348, ±0.522 GeV 2 are calculated. These correspond to k d = ±1, ±2, ±3, ±4, ±6 in Eq. 5. The zeromomentum projected correlation functions contain spinup and spin-down components.
It is important to note that these configurations are electro-quenched; the field exists only for the valence quarks of the hadron. To include the background field at configuration generation time is possible [20] but requires a separate Monte Carlo simulation for each field strength and is hence prohibitively expensive. Separate calculations also destroy the advantageous correlations between the field strengths used when constructing the ratio in Eq. 15. An alternative is to use a reweighting procedure on the gauge field configurations [21] for the different field strengths B, but this is not performed here.

V. FITTING
The energy shift at each field strength has the form specified by Eq. (12) and as such we fit with a quadratic term, Figures 6 and 7 show fits for the neutron energy shift with a smeared source and a U (1) Landau mode sink projection. For the first time, clear plateaus are present in this difficult to obtain quantity. It is required that a plateau be present at each of the three non-zero field strengths in order to proceed to the next stage. The plateau does not occur until t = 24, this region is a common starting point across the heavier masses. The primary cause for this late plateau onset time is the zero-field correlator which has fundamentally different physics. As such its potential excited state behaviour is different to that of the background field correlators. Plateaus only form once both correlators have decayed to the ground state.
The fit performed is as a function of k d , the integer magnetic flux quanta in Eq. (5), Here c 2 is the fit parameter, and has units of GeV as these are the units of δE(k d ). In order to convert this fit parameter to the physical units of magnetic polarisability, fm 3 , Eq. (5) is used to produce the transformation Here α is the fine structure constant, α ≈ 1/137. The quadratic fitting process uses only energy shifts which have the same spatial modulation of Eq. (10). As was seen in Figures 2 and 3, the spatial extent affects the coupling to the energy eigenstates. It is hence important that when we fit, we use the optimal sink projection. This is achieved using a simultaneous investigation of the spatial extent, σ z and field strength. For a specified spatial modulation to be suitable it must provide early isolation of the eigenstate at each field strength. This isolation is visible in the long plateaus of Figures 2 and 3. This is already a strong constraint on the sink choice but in order to determine where to fit energy shifts for the quadratic fit in B a further constraint is needed. This constraint comes from considering the constant plateau fits to the energy shift at all field strengths. By considering all possible fit windows, we select fit windows where  good plateau behaviour exists for all field strengths simultaneously. Good plateau behaviour is characterised by a fit χ 2 dof of less than 1.2. This proecess of requiring good plateau behaviour at each field strength simultaneously dramatically reduces the number of possible fit windows. In particular, it is often difficult to obtain acceptable energy shift plateaus for the largest field strength considered.
The final constraint on the fitting process comes from the quadratic fit itself. This fit must also be acceptable having a χ 2 dof ≤ 1.2. If multiple fit windows remain after this process, the one with the longer time extent and χ 2 dof 's closest to one are preferred. Once the specific quadratic fit has been chosen, the magnetic polarisability, β is extracted from the quadratic coefficient of the fit. In order to test the presence of higher order terms in the energy shift of Eq. (17), a quartic term is also considered. It is found that the quartic term is not needed in order to fit the energy shifts well with acceptable χ 2 dof . Using the sink eigenmode-projection technique at each quark mass, it is possible to extract magnetic polarisabilities from the fits to the constant energy shift plateaus as a function field strength. Results are presented in Table  I for the magnetic polarisability of the neutron at each quark mass.

VI. CHIRAL EXTRAPOLATION
Chiral effective field theory (χEF T ) is an important tool for connecting lattice results to the physical point. The extrapolation here follows that in Ref. [22] and is summarised briefly below.
The lattice QCD results do not incorporate contributions from photons coupling to sea-quark loops which form the meson dressing of χEF T -they are electroquenched. Thus it is necessary to model the corrections associated with these effects, which is done using partially quenched χEF T for the neutron.
Finite volume effects are considered by replacing the continuum integrals of the chiral expansion with sums over the momenta available on the lattice. We note that the lattice volume varies slightly across the four lattice data points available due to our use of the Sommer scale.
The leading-order loop contributions are calculated using the heavy-baryon approximation following [23,24]. We consider the chiral expansion where the residual series coefficients, a 0 , a 2 , are constrained by our lattice QCD results, following corrections to infinite volume. When the finite-range regularisation approach is used, the leading non-analytic contribution to the chiral expansion comes entirely from the diagram shown in Figure 9 for the above parameterisation of β n . Figure 10 shows transitions to a Delta baryon at leading-order; for a finite mass-splitting ∆ this contributes a log term rather than 1/m π [24].
The three-dimensional integral forms are [22] β (π N ) m 2 π = e 2 4 π 1 288 π 3 f 2 π χ N d 3 k k 2 u 2 (k, Λ) Here ω k = k 2 + m 2 π is the energy carried by the pion which has three-momentum k, ∆ is the mass splitting between the Delta baryon and the nucleon, ∆ ≡ M ∆ − M N = 292 MeV and the pion decay constant is taken as f π = 92.4 MeV. The standard coefficients for full QCD are modified to account for partial quenching effects as explained in Ref. [22]. The coupling constants take the values g A = 1.267 and C = −1.52.
The dipole regulator, u (k, Λ) , ensures that only soft momenta flow through the effective field theory degrees of freedom. The value Λ = 0.80 GeV is adopted [25][26][27][28][29] as this regulator mass defines a pion cloud contribution to masses [26], magnetic moments [27] and charge radii [25], enabling a correction to the pion cloud encountered in unquenching. This particular choice of regulator mass defines a neutron core contribution which does not differ significantly between partially quenched QCD and full 2+1 flavour QCD. This connection relies on the assumption that the core is insensitive to sea-quark loop contributions [30].

A. Results
We proceed by calculating the integrals of Eqs. (21) and (22) in finite and infinite volume and using the differences at each quark mass considered on the lattice to correct the finite volume lattice-QCD results to infinite volume. In addition, the chiral coefficients of Eqs. (23) and (24) are varied. In the finite volume, particallyquenched coefficients, omitting the contributions of direct sea-quark-loop photon interaction contributions, are considered [22]. In infinite volume, the full-QCD coefficients of Eqs. (23) and (24)   volume and sea-quark-loop contribution corrections are incorporated. These corrections are illustrated in Figure  11 by the (orange) "Full-QCD Infinite-Volume Results" next to the original (blue) "Lattice Results". With the lattice results corrected to infinite-volume full-QCD results, the fit parameters a 0 and a 2 are constrained by fitting Eq. (20) to the corrected lattice results. Once determined, any volume can be considered. Figure 11 also shows chiral extrapolations for a range of volumes. Large box sizes are requried in order to obtain an extrapolation close to the infinite volume value at the physical point. These extrapolations have also been corrected for electro-quenching effects and hence are a guide to future lattice QCD simulations which include the background field effects on the disconnected sea-quarks.
Following the procedure summarised above [22], the physical polarisability can be obtained from Eq. (20) with m π = m phys π = 140 MeV. We find β n = 2.02(16)(8) × 10 −4 fm 3 at the physical point. The stated uncertainties are derived from the statistical error of the fit parameters of the lattice QCD results and the systematic uncertainty from variation of the regulator Λ over the range 0.7 ≤ Λ ≤ 0.9 GeV respectively. A comparision between this result and the experimental data is provided in Figure 12. Our calculation is in good agreement with a number of the experimental results and poses an interesting challenge for greater experimental precision. Similarly, progress in experimental measurement would drive further lattice QCD and chiral effective field theory work.
These lattice results use a single lattice-spacing and as such it is not possible to quantify an uncertainty associated with taking the continuum limit. However as a non-perturbatively improved clover fermion action is used the O a 2 corrections are expected to be small relative to the uncertainties already presented. It is anticipated that there is some degree of additive quark mass renormalisation due to the interaction of the background field with the Wilson term in the fermion action [32] but the extent to which this effect remains with the clover fermion action is under investigation [33].

VII. CONCLUSION
The neutron magnetic polarisability has been calculated using a novel approach in which asymmetric operators are used at the source and sink. The use of gaugeinvariant Gaussian smearing at the source encapsulates the dominant QCD dynamics, while a gauge-fixed U (1) two-dimensional eigenmode projection technique is used at the sink to encode the Landau level physics resulting from the presence of the uniform magentic field. A systematic exploration of the parameter space was used to optimise operators that couple efficiently to the neutron ground state in a magnetic field. The use of this Landau mode projection at the sink has for the first time enabled the fitting of plateaus in the magnetic polarisability energy shifts.
Calculations at several pion masses have enabled the use of heavy baryon chiral effective field theory to relate the lattice QCD results to experiment. This enables us to make a theoretical prediction for the neutron magnetic polarisability of β n = 2.02(16)(8) × 10 −4 fm 3 . This prediction is founded by ab initio lattice QCD simulations, with effective field theory techniques used to account for the finite volume of the lattice, disconnected sea-quarkloop contributions, and extrapolate to the light quark masses of Nature. The resulting value is in agreement with the current experimental estimates and presents an interesting challenge for greater experimental precisioin.
A range of finite-volume extrapolations are performed in the framework of chiral effective field theory. As these curves incorporate the contributions of sea-quark loops they are presented as a guide to future lattice QCD simulations. At the physical pion mass, the 7 fm curve still differs from the infinite-volume prediction by 6%.
Our result does not directly incorporate the sea-quark loop effects of the magnetic field in the lattice simulation, which would require a separate Monte Carlo ensemble for each value of B considered and as such is prohibitively expensive. We also note that in this study we have not considered the effects from the B-dependent additive quark mass renormalisation that arises due to the use of Wilson-type fermions [32]. The detailed impact of this effect on clover fermions has not been studied and will be addressed in a separate work that is in preparation [33]. While the effect can be observed in pion correlators, our preliminary tests indicate it is hidden by the large statistical fluctuations inherent in baryon correlators.
It is interesting to consider that, in obtaining the magnetic polarisability, we want to work with small B field strengths in order to make use of the perturbative energy expansion for the neutron. This means we are the confining phase of QCD, such that quarks cannot have individual Landau levels. Nonetheless, the success of our Landau mode sinks indicates that the effects of the magnetic field on the quark distribution in the neutron are significant.
Future work will examine the effect of using a gauge-covariant sink projection based on the eigenmodes of the two-dimensional U (1)×SU (3) Laplacian [15]. This alters the Landau level structure, breaking the degeneracy and mixing different levels. Indeed, a recent finite temperature study using the staggered quark formulation finds that the contribution of the lowest Landau level eigenmodes remains important even after QCD interactions are introduced [14], further motivating our investigation of the effectiveness of an eigenmode projected sink which is aware of the QCD gauge field.