Jean-Gabriel Minonzio
email: jean-gabriel.minonzio@uv.cl

Chao Han

Didier Cassereau

Quentin Grimal

In vivo pulse-echo measurement of apparent broadband attenuation and Q factor in cortical bone: A preliminary study

Keywords: cortical bone, broadband ultrasonic attenuation, Quantitative Ultrasound, in vivo measurement, Orthogonal Matching Pursuit, sparse reconstruction, Q factor

Quantitative UltraSound (QUS) methods have been introduced to assess cortical bone health at the radius and tibia through the assessment of Cortical Thickness (Ct.Th), Cortical Porosity (Ct.Po) and bulk wave velocities. Ultrasonic attenuation is another QUS parameter which is not currently used. We assess the feasibility of in vivo measurement of ultrasonic attenuation in cortical bone with a broadband transducer with 3.5 MHz-center frequency. Echoes from the periosteal and endosteal interfaces were fitted with Gaussian pulses using sparse signal processing. Then, the slope of the Broadband Ultrasonic Attenuation (Ct.nBUA) in cortical bone and quality factor Q -1

were calculated with a parametric approach based on the center-frequency shift. Five human subjects were measured at the one-third distal radius with pulse-echo ultrasound, and reference data was obtained with high-resolution X-ray peripheral computed tomography (Ct.Th and Cortical volumetric Bone Mineral Density, Ct.vBMD). Ct.Th was used in the calculation of Ct.nBUA while Q -1 11 is obtained solely from ultrasound data. The values of Ct.nBUA

Introduction

Osteoporosis fracture risk is currently assessed using Dual energy Xray Absorptiometry (DXA) in order to assess areal Bone Mineral Density (aBMD). However, DXA has strong limitations, in particular it lacks sensitivity [START_REF] Briot | FRAX R : Prediction of Major Osteoporotic Fractures in Women from the General Population: The OPUS Study[END_REF][START_REF] Siris | Absolute Versus Relative Fracture Risk[END_REF] and is not appropriate to monitor cortical bone [START_REF] Choksi | The challenges of diagnosing osteoporosis and the limitations of currently available tools[END_REF]. Cortical bone, the dense tissue that forms the outer shells of the bones, represents about 80% of the human skeleton mass and plays an important role in the skeletal mechanical stability [START_REF] Holzer | Hip Fractures and the Contribution of Cortical Versus Trabecular Bone to Femoral Neck Strength[END_REF][START_REF] Zebaze | Intracortical remodelling and porosity in the distal radius and post-mortem femurs of women: a cross-sectional study[END_REF][START_REF] Bala | Cortical Porosity Identifies Women With Osteopenia at Increased Risk for Forearm Fractures[END_REF]. Aging and bone pathologies are associated with cortical thinning [START_REF] Nishiyama | Postmenopausal Women With Osteopenia Have Higher Cortical Porosity and Thinner Cortices at the Distal Radius and Tibia Than Women With Normal aBMD: An In VivoHR-pQCT Study[END_REF] and weakening of the bone material mechanical quality reflected in an increase of Cortical Porosity (Ct.Po) [START_REF] Kral | Increased cortical porosity and reduced cortical thickness of the proximal femur are associated with nonvertebral fracture independent of fracture risk assessment tool and garvan estimates in postmenopausal women[END_REF] or a decrease of Cortical volumetric Bone Mineral Density (Ct.vBMD (CT.vBMD) [START_REF] Ostertag | Multiscale and multimodality computed tomography for cortical bone analysis[END_REF][START_REF] Neto | Vitamin d deficiency is associated with cortical bone loss and fractures in the elderly[END_REF].

Several quantitative ultrasound (QUS) approaches have been introduced

to assess cortical bone health at the radius and tibia. Some aim at assessing Cortical Thickness (Ct.Th) assuming a nominal value of ultrasound velocity using pulse-echo [START_REF] Karjalainen | Ultrasonic assessment of cortical bone thickness in vitro and in vivo[END_REF], axial transmission [START_REF] Moilanen | Ultrasonic guided waves in bone[END_REF], or through transmission measurements [START_REF] Sai | Novel ultrasonic bone densitometry based on two longitudinal waves: significant correlation with pQCT measurement values and age-related changes in trabecular bone density, cortical thickness, and elastic modulus of trabecular bone in a normal Japanese population[END_REF]. Other approaches are designed for a combined estimation of Ct.Th and Ct.Po or bulk wave velocities using axial transmission measurements of several guided wave modes [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: A pioneering study on in vitro cortical bone samples[END_REF][START_REF] Minonzio | Ultrasound-Based Estimates of Cortical Bone Thickness and Porosity Are Associated With Nontraumatic Fractures in Postmenopausal Women: A Pilot Study[END_REF] or adaptative pulse-echo imaging with a transducer array [START_REF] Renaud | In vivo ultrasound imaging of the bone cortex[END_REF][START_REF] Renaud | Measuring anisotropy of elastic wave velocity with ultrasound imaging and an autofocus methodapplication to cortical bone[END_REF]. Ultrasonic attenuation is another QUS parameter which has been exploited for several decades in soft tissues [START_REF] Mamou | Quantitative Ultrasound in Soft Tissues[END_REF] and trabecular bone [START_REF] Langton | The measurement of broadband ultrasonic attenuation in cancellous bone-a review of the science and technology[END_REF]; however, until now, it has received little consideration in cortical bone QUS.

Ex vivo studies on cuboid specimens have established that ultrasonic attenuation of bulk waves in cortical bone is related to mass density [START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia middiaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF] and to Ct.vBMD [START_REF] Sasso | Dependence of ultrasonic attenuation on bone mass and microstructure in bovine cortical bone[END_REF]. Also, simulations have suggested that the scattering of ultrasound by the cavities of the pore network is one important mechanism of attenuation [START_REF] Yousefian | The effect of pore size and density on ultrasonic attenuation in porous structures with mono-disperse random pore distribution: A two-dimensional in-silico study[END_REF][START_REF] Yousefian | Frequency-dependent analysis of ultrasound apparent absorption coefficient in multiple scattering porous media: application to cortical bone[END_REF][START_REF] Iori | Estimation of Cortical Bone Microstructure from Ultrasound Backscatter[END_REF]. It follows that, in addition to ultrasonic velocities measured with QUS approaches [START_REF] Grimal | Quantitative ultrasound assessment of cortical bone properties beyond bone mineral density[END_REF], attenuation could be indicative of the mechanical quality of cortical bone as it is related to Ct.Po and Ct.vBMD.

There has been a few attempts to measure attenuation in ex vivo bone specimens. [START_REF] Zheng | Spectral ratio method to estimate broadband ultrasound attenuation of cortical bones in vitrousing multiple reflections[END_REF] measured in pulse echo mode a bovine femur specimen and estimated the slope of the frequency-dependent attenuation coefficient (also refered to as "spectral ratio method") to estimate the so-called Cortical normalized Broadband Ultrasonic Attenuation (Ct.nBUA). They later used a parametric approach introduced by [START_REF] Kuc | Parametric Estimation of the Acoustic Attenuation Coefficient Slope for Soft Tissue[END_REF] which related attenuation to the shift of the center frequency (also refered to as "peak frequency method") [START_REF] Zheng | Broadband ultrasound attenuation measurement of long bone using peak frequency of the echoes[END_REF]. [START_REF] Dencks | Modelbased estimation of quantitative ultrasound variables at the proximal femur[END_REF] also used this parametric approach and measured nBUA in proximal femurs in throughtransmission. However, these values of nBUA can hardly be interpreted in terms of bulk cortical bone material properties because ultrasound propagated along a complex path through both cortical and trabecular bone. As far as we know, in vivo measurements of attenuation in cortical bone have not yet been reported.

The aim of this paper was to assess the feasibility of in vivo measurement of attenuation in cortical bone. We have conducted a preliminary study on the radius of five human subjects. Ultrasound echoes stemming from the normal-incidence reflection on the outer (periosteal) and inner (endosteal) cortical bone interfaces were recorded and processed with Orthogonal Matching Pursuit (OMP) to retrieve the time delay, the temporal echo width and the frequency shift of the center frequency from which Ct.nBUA and the quality factor Q were calculated with a parametric approach [START_REF] Kuc | Parametric Estimation of the Acoustic Attenuation Coefficient Slope for Soft Tissue[END_REF]. The quality factor is introduced as a quantity related to the dissipation of energy which can be measured without the need of the knowledge of the bulk wave velocity. Ultrasound parameters Ct.nBUA and Q were compared to reference values of Ct.vBMD of each subject which were obtained with High-Resolution X-ray peripheral Quantitative Computed Tomography (HR-pQCT). The results are of interest for the development of future multimodal cortical bone QUS approaches estimating bone structural (thickness)

and material (porosity, velocities, attenuation) properties for the evaluation of bone health.

Method

Extraction of echoes

The signal received in the pulse-echo ultrasound measurement is modeled as a sum of two Gaussian pulses s i (t), or Gabor functions, i.e., the product of a Gaussian function with a complex sinusoid [START_REF] Demirli | Model-based estimation of ultrasonic echoes. Part II: Nondestructive evaluation applications[END_REF])

y(t) = 2 i=1 s i (t) + n(t) = 2 i=1 A i exp[- (t -t i ) 2 2σ 2 t i ] cos[2πf i (t -t i )] + n(t), (1)
where σ t i is the standard deviation of the temporal Gaussian function, f i is the central frequency, t i is the group delay and A i is the amplitude of i-th echo. The error between the model and the measured signal, including noise, is n(t). The first and second echoes correspond to the reflections on the outer (periosteal) and inner (endosteal) bone surfaces, respectively.

For further use, we write the temporal Fourier transform of one echo as

S i (f ) = A i √ 2πσ f i exp[- (f -f i ) 2 2σ 2 f i ] exp[-j2πf t i ], (2)
where the standard deviation of the Gaussian function σ f i satisfies 2πσ f i = 1/σ t i . The -6dB frequency bandwidth, or full width at half maximum, is related to σ f i , and is equal to 2 2 ln(2)

σ f i ≈ 2.35σ f i .
The echoes are isolated using a sparse signal processing method in the time domain which provides the quantities σ t i , f i , t i , and A i (i = 1, 2). The method is detailed in Appendix ; shortly, the signal model, given by (1), discretized in time can be represented as

y = Dx + n, (3)
where y and n are N t × 1 vectors corresponding to the sampling of y(t) and n(t) at N t discrete time points, respectively. Likewise, x is a N m × 1 vector collecting the N m relative amplitudes of the echoes. In a sparse point of view, only a few elements of this vector are non zero. Finally, D is a N t × N m matrix corresponding the so-called dictionary. Each column

d(t; σ tm , f m , t m ), with m = 1 • • • M , of D is one of the M Gabor functions discretized at
N t time points, defined by the set of parameters (σ tm , f m , t m ). For a given measurement y(t), the problem amounts to determine the two non-zero components of x, which indices yield the set of parameters corresponding to the two echoes. This is done by sparse reconstruction with Orthogonal Matching Pursuit (OMP) [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF]. The details of the construction of the dictionary and of the implementation of OMP can be found in Appendix.

Measurement of attenuation: Theory

The cortical bone layer is modeled locally as a plate of thickness Ct.Th.

We assume that the temporal signals can be decomposed into monochromatic plane waves propagating in bone at normal incidence on the plate surfaces.

The complex wavenumber is denoted (k + iα), where α is the imaginary part and is assumed to present a linear attenuation with frequency in cortical bone, i.e., α = β f [START_REF] Minonzio | Impact of attenuation on guided mode wavenumber measurement in axial transmission on bone mimicking plates[END_REF]. This linear approximation, or first order Taylor expansion, is valid around a central frequency f 0 as long as the frequency deviation ∆f is narrow, i.e., ∆f /f 0 << 1, even if the attenuation variation for larger frequency is not linear [START_REF] Szabo | Causal theories and data for acoustic attenuation obeying a frequency power law[END_REF][START_REF] Yousefian | Frequency-dependent analysis of ultrasound apparent absorption coefficient in multiple scattering porous media: application to cortical bone[END_REF]. Note that the real part k of the wavenumber is assumed to be frequency independent within the considered bandwidth. With this plane wave model, the modulus of the ratio of the spectra of the two echoes (from the periosteal and endosteal surfaces) is proportional to exp(-βf 2Ct.T h) [START_REF] Zheng | Broadband ultrasound attenuation measurement of long bone using peak frequency of the echoes[END_REF]. Note that the ultrasound signal does not need to be corrected for the overlying soft tissues, as we are studying the ratio between endosteal and periosteal echoes.Both echoes are indeed equally affected by the propagation within the soft tissue layer. Thus, the ratio only depends on the propagation inside the cortical bone layer.

Following [START_REF] Kuc | Parametric Estimation of the Acoustic Attenuation Coefficient Slope for Soft Tissue[END_REF], the attenuation coefficient β can be estimated from pulse-echo measurements with a parametric approach from the shift of the center frequency of the Gaussian pulse. Assuming that the central frequency variation ∆f = f 1 -f 2 is small compared to the central frequency (i.e., σ f remains unchanged), the coefficient β writes [START_REF] Kuc | Parametric Estimation of the Acoustic Attenuation Coefficient Slope for Soft Tissue[END_REF][START_REF] Narayana | On the validity of the linear approximation in the parametric measurement of attenuation in tissues[END_REF])

β = 1 2Ct.T h ∆f σ 2 f , (4)
with the frequency in MHz and distance in mm, β is in Np.mm -1 .MHz -1 .

This parametric estimation is well adapted in both transmission [START_REF] Kuc | Parametric Estimation of the Acoustic Attenuation Coefficient Slope for Soft Tissue[END_REF] and reflection [START_REF] Kuc | Estimating acoustic attenuation from reflected ultrasound signals: Comparison of spectral-shift and spectral-difference approaches[END_REF] and has been successfully applied to bone on ex vivo specimens in both configurations [START_REF] Dencks | Modelbased estimation of quantitative ultrasound variables at the proximal femur[END_REF][START_REF] Zheng | Spectral ratio method to estimate broadband ultrasound attenuation of cortical bones in vitrousing multiple reflections[END_REF][START_REF] Zheng | Broadband ultrasound attenuation measurement of long bone using peak frequency of the echoes[END_REF]. Finally, the cortical broadband ultrasound attenuation in dB.cm -1 .MHz -1 , is obtained as

Ct.nBUA = 10 20 ln(10) β ≈ 86.9β.

(5)

The above equations indicate that the experimental determination of nBUA or β requires the knowledge of Ct.Th. This can be obtained from the X-ray computed tomography scan of the bone, or, alternatively, from ultrasound signals as

Ct.Th = v 11 ∆t 2 , (6)
where ∆t = t 2 -t 1 , providing the longitudinal bulk wave velocity v 11 (in mm.µs -1 ) is known. The index "11" refers to the radial bone direction commonly denoted direction 1 in previous studies [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: A pioneering study on in vitro cortical bone samples[END_REF][START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia middiaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF]. In a first approach, a nominal value of v 11 may be assumed, i.e., the velocity is supposed to be known and identical for all subjects [START_REF] Karjalainen | Ultrasonic assessment of cortical bone thickness in vitro and in vivo[END_REF][START_REF] Grimal | Quantitative ultrasound assessment of cortical bone properties beyond bone mineral density[END_REF]. Thus, the β coefficient may be rewritten without reference to Ct.Th, using equations ( 4) and ( 6), as

β = 2 π v 11 σ t ∆t ∆f σ f . (7)
The quality factor Q is a dimensionless parameter, classically used to describe the resonance of a resonator, usually defined by the ratio between the central frequency and the frequency bandwidth. High Q values correspond to low attenuation. It is possible to define a quality factor related bulk wave velocity v 11 as Q -1 11 = (C 11 ) (C 11 ) [START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia middiaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF], where C 11 is the complex elastic coefficient related to v 11 . Interestingly, as we show below, Q -1 11 does not require the knowledge of the bulk wave velocity. In case of weak attenuation, i.e., Q -1 11 1 or α k, which is usually satisfied in cortical bone at low frequency, i.e., less than a few MHz, [START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia middiaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF], the quality factors writes

Q -1 11 ≈ 2α k ≈ βv 11 π . (8)
Which can advantageously be rewritten

Q -1 11 ≈ 1 πσ 2 f ∆f ∆t , (9)
which expression does not depend on the bulk wave velocity v 11 but only on the parameters ∆t, ∆f and σ f which can be extracted from the measured signals. Thus, evaluating Q -1 11 to measure attenuation could be an advantage as this does not require the knowledge of the bulk wave velocity.

In vivo ultrasound measurements

This study has been approved by the ethical committee of the Committees Those spectra were obtained by the temporal Fourier transforms of segments s 1 (t) and s 2 (t), indicated with thick lines, corresponding to 1.2 µs from each side of the envelop maxima [START_REF] Karjalainen | Ultrasonic assessment of cortical bone thickness in vitro and in vivo[END_REF].

Finally, on Figure 1(c), one can observe that the variation of the spectrum ratio ln(|S 2 (f )/S 1 (f )|) is in agreement with the slope ∆f /σ 2 f of the linear approximation [START_REF] Kuc | Parametric Estimation of the Acoustic Attenuation Coefficient Slope for Soft Tissue[END_REF]. Note that the spectral ratio method is based on the slope evaluation, while the peak frequency method is based on Kuc's formula. Processing of the signal with OMP yields the parameters (σ f 1 , f 1 , t 1 ) and (σ f 2 , f 2 , t 2 ) of the echoes which were used to calculate Ct.nBUA (eq. 5) using the cortical thickness value obtained from HR-pQCT (see section 1.4) (eq. 4) and Q -1 11 (eq. 9). The extraction of the parameters (σ t 1 , f 1 , t 1 ) and (σ t 2 , f 2 , t 2 ) is illustrated on Figure 1

HR-pQCT

The left radius of each subject was scanned with high-resolution X-ray peripheral computed tomography (HR-pQCT, XtremeCT, Scanco Medical, Brüttisellen, Switzerland) using the standard manufacturer in vivo acquisition parameters (60 kVp, 1000 mA, 100 ms integration time, voxel size 82 µm). The measurements were approximately performed at the one-third distal radius site. To ensure that HR-pQCT imaging is site-matched with ultrasound measurement, we used an initial large view image allowing to center the scanning windows at 7 cm from the radial styloid corresponding to the ultrasound measurement point. The acquired volume corresponded to 9.02 mm along the bone axis (110 cross-sectional slices). The central slice is illustrated in Figure 3 for each subject. The region of interest (ROI), site-matched with the ultrasound measurement site, is highlighted. Methods used to process the CT data have been previously described in detail [START_REF] Laib | In vivo high resolution 3D-QCT of the human forearm[END_REF][START_REF] Boutroy | In Vivo Assessment of Trabecular Bone Microarchitecture by High-Resolution Peripheral Quantitative Computed Tomography[END_REF][START_REF] Macneil | Accuracy of high-resolution peripheral quantitative computed tomography for measurement of bone quality[END_REF][START_REF] Ostertag | Multiscale and multimodality computed tomography for cortical bone analysis[END_REF]. Briefly, cortical limits were determined using a threshold-based The relationships between ultrasound parameters and reference parameters obtained with HR-pQCT were assessed using linear correlations. We report Pearson's correlation coefficient and linear regression equations.

Results

The mean and standard deviation, calculated over the ten retained measurements, of the estimated parameters ∆t, ∆f , σ f for each subject are reported in Table 1. The next last column corresponds to Q -1 11 obtained combining previous parameters using equation (9). Finally, values of Ct.Th, ∆f and σ f were used to compute Ct.nBUA according to equations ( 4) and

(5) (last column of Table 1). HR-pQCT parameters Ct.Th and Ct.vBMD are reported in Table 2.

High correlations were observed between Q -1 11 and Ct.vBMD, obtained by HR-pQCT (R 2 = 0.93, p < 0.01, RMSE= 0.76%, Fig. 4) and between

Ct.nBUA and Ct.vBMD (R 2 = 0.92, p < 0.01, RMSE = 0.56 dB.MHz -1 .cm -1 , Fig. 5). The linear regression equation was Ct.nBUA = -34.1 × Ct.vBMD + 45.6, where Ct.nBUA is expressed in dB.MHz -1 .cm -1 and Ct.vBMD in g.cm -3 . Likewise, the second linear regression equation was Q -1 11 = -49.2 × Ct.vBMD + 64.9, where Q -1 11 is expressed in % and Ct.vBMD in g.cm -3 . The corresponding longitudinal bulk wave velocities v 11 , ranging from 3.2 to 3.9 mm.µs -1 , are given in Table 3 for comparison with other studies.

Discussion

This work considered pulse-echo measurements to measure the attenuation of signals from ultrasonic waves propagating in cortical bone at the radius in vivo. The two echoes stemming from normal-incidence reflections on the periosteal and endosteal bone surfaces were modeled as the sum of two elementary waveforms (Gabor functions), and the parameters of these waveforms (time delay, central frequency and frequency bandwidth) were recovered using sparse signal processing (OMP). Attenuation was assessed in two ways: (1) as Ct.nBUA with a parametric method using OMP parameters combined with Ct.Th from HR-pQCT ; and (2) as a quality factor Q -1 11 calculated from OMP parameters only.

The values of Ct.nBUA (6.7±2.2 dB.MHz. -1 .cm -1 ) are in the range of the attenuation values reported in the literature on bovine and human femur (Table 4). To which extent Ct.nBUA and Q -1 11 reflect the true (intrinsic) ultrasonic attenuation in cortical bone remains to be investigated. Such quantity can only be evaluated ex vivo using dedicated experimental conditions in order to minimize the effect of diffraction and other losses unrelated to attenuation within the tissue.

We found that Ct.nBUA was strongly correlated to volumetric bone mineral density (Ct.vBMD) measured from X-ray attenuation (R 2 = 0.92, Fig. 5).

This finding is in line with the ex vivo results of Sasso et al. [START_REF] Sasso | Dependence of ultrasonic attenuation on bone mass and microstructure in bovine cortical bone[END_REF] who reported a correlation between Ct.nBUA and Ct.vBMD in bovine bone measured ex vivo (R 2 = 0.57, p < 10 -5 , RMSE=1.6). These authors also reported the linear fit equation between Ct.nBUA and Ct.vBMD Quality factor Q -1 11 (Table 1) values associated to the propagation of a longitudinal wave can be compared to shear mode quality factor Q -1 44 measured from the first resonance peak (falling in the range [100-300 kHz]) of a cuboid bone specimen [START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia middiaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF]. In the latter study, the average Q -1 44 was 3.5% to be compared to a mean value of 8.6% in the present study. These values compare well although the comparison should be made with caution because (i) Q -1 44 was obtained at a much lower frequency (one order of magnitude); (ii) the polarization is different and a larger attenuation is expected for shear waves; (iii) the collection of samples used in [START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia middiaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF] includes low density (high porosity) samples which may not be representative of the bone of healthy volunteers in the present study.

Ct.nBUA was calculated from measured ultrasound parameters and a value of Ct.Th obtained from the HR-pQCT image of each subject. If only ultrasound data is available, the quality factor Q -1 11 can be calculated as it does not depend on the bulk wave velocity (Eq. 9). Our results suggest that Q -1 11 could be of clinical interest as it is correlated with Ct.vBMD (R 2 = 0.92, Fig. 4). Furthermore, [START_REF] Fan | Cortical bone viscoelastic damping assessed with resonant ultrasound spectroscopy reflects porosity and mineral content[END_REF] have shown that the shear wave quality factor is correlated to porosity (R 2 =0.53), and it is reasonable to infer that

Q -1
11 measured in this study is also related to porosity.

Signal parameters were estimated in the ultrasound signal with OMP, which has several advantages. The segmentation of the two echoes in the time domain was performed automatically in a robust manner, which can be an advantage in case of overlapping of echoes (which can occur for thin cortices). Also, because the Gabor function offers a reliable parametrization of the echoes, the calculation of Fourier domain parameters is done avoiding a Fourier transform of the signal which can be polluted by the choice of the time window for time segmentation.

The parametric method used to estimate attenuation relies on the measurement of ∆f . In equation ( 7), β is written in terms of temporal and frequency shifts divided by their associated standard deviations, leading to normalized shifts or Z-scores. Using equations ( 8), Q 11 -1 can also be written in term of these ratios as ∆f

σ f = 1 2 Q 11 -1 ∆t σt .
As the ratio ∆t/σ t is of the order of 1 and Q -1 11 is small compared to 1, one can expect the ratio ∆f /σ f to be also small compared to 1 in accordance with the weak attenuation hypothesis.

Thus, special attention should be paid to the evaluation of ∆f . Indeed, it can be observed in Table 1 that the largest relative uncertainties, up to 70%, are obtained for the frequency ∆f parameter, while the uncertainties on the temporal ∆t parameter are about a few percents. This observation should be taken into account in order to propose a robust clinical measurement.

The calculation of Ct.nBUA requires the knowledge of cortical thickness Ct.Th, which was obtained from HR-pQCT in the present proof-of-concept study. Our method to measure attenuation could however be implemented together with other ultrasound sequences and signal processing providing speed of sound and cortical thickness as described in [START_REF] Renaud | In vivo ultrasound imaging of the bone cortex[END_REF] or

Nguyen [START_REF] Nguyen Minh | Estimation of Thickness and Speed of Sound in Cortical Bone Using Multifocus Pulse-Echo Ultrasound[END_REF].

This study has a number of limitations. (1) Only 5 healthy subjects were considered. The range of cortical bone properties (thickness, material properties) may not be representative of the general population or subjects with bone pathologies.

(2) The endosteal interface of cortical bone of young healthy subjects is known to be quasi-plane and regular, whereas it is likely to be irregular and discontinuous for elderly subjects or patients with bone diseases as a result of age-or disease-related bone deterioration and trabecularization [START_REF] Zebaze | Intracortical remodelling and porosity in the distal radius and post-mortem femurs of women: a cross-sectional study[END_REF]. Such irregular interface may give rise to a wave pattern more complex than the specular echo observed in the young subjects in the present study. Also, as the attenuation is expected to increase with porosity, the signal from the endosteal interface in osteoporostic or old patients should be weaker compared to the signals processed in this study.

The applicability of our method to assess aged or diseased bone remains to be assessed. (3) Because the probe was handheld, some pulse echo acquisitions had to be discarded. We used a quality criterion provided by OMP processing based on the value of the scalar product between the recorded signal and the dictionary. This criterion is believed to be robust as it relies on the expected shape of the echo waveform. Nevertheless, reproducibility of the method should be assessed in future studies on a larger number of subjects. Moreover, all measurements have been carried out at the same central frequency of 3.5 MHz. This frequency may need to be adapted for smaller or larger thicknesses in order to remain with an endosteal echo clearly separated and not too much attenuated. As attenuation in cortical bone may have a nonlinear dependence on frequency [START_REF] Yousefian | Frequency-dependent analysis of ultrasound apparent absorption coefficient in multiple scattering porous media: application to cortical bone[END_REF], attenuation values obtained at different frequencies may not be directly compared with the values obtained in the present study.

The results of this study suggest that ultrasonic quantities (Ct.nBUA and [START_REF] Mohanty | Artificial neural network to estimate micro-architectural properties of cortical bone using ultrasonic attenuation: A 2-d numerical study[END_REF], machine learning [START_REF] Minonzio | Automatic classifying of patients with non-traumatic fractures based on ultrasonic guided wave spectrum image using a dynamic support vector machine[END_REF] or deep learning [START_REF] Li | Deep learning analysis of ultrasonic guided waves for cortical bone characterization[END_REF] will be investigated in order to improve the robustness of the approach.

Q -1 11 ) related
d(t, Θ) = ζ Θ exp[-s(t -τ ) 2 ] exp[j2πf (t -τ )] (.1) Θ = [s, τ, f ] ( .
M = L × S × K, with M N t .
The dictionary can be represented as

D =      d(t 1 , Θ 1 ) d(t 1 , Θ 2 ) • • • d(t 1 , Θ M ) d(t 2 , Θ 1 ) d(t 2 , Θ 2 ) • • • d(t 2 , Θ M ) . . . . . . . . . . . . d(t Nt , Θ 1 ) d(t Nt , Θ 2 ) • • • d(t Nt , Θ M )      ( .3)
where

Θ m = [s l , τ s , f k ], l ∈ [1, L], s ∈ [1, S], k ∈ [1, K], m ∈ [1, M ].
The ranges of variations of s, τ and k are defined as follows. The range of frequency f can be obtained from the bandwidth of the received signal, as discussed in [START_REF] Mor | A matching pursuit method for approximating overlapping ultrasonic echoes[END_REF]. Similarly, the range for the bandwidth factor s can be deduced from the bandwidth of the emitted signal. The time delays τ of the different echoes are searched around the peaks of the envelope of the received signal. Overall, the selection of the parameter intervals is quite flexible due to the robustness and efficiency of the OMP method. Using the dictionary, the signal model (1) can be represented as

y = Dx + n (.4)
where y and n are the sampling of y(t) and n(t) at discrete time points respectively, x is the amplitude vector corresponding to each column in D.

This dictionary can be interpreted as an extension of the classical Fourier transform. In the case of the temporal Fourier transform, Θ is reduced to the frequency f and a element of the dictionary matrix writes as d(t, Θ) = d(t, f ) = exp(-j2πf t) and previous equation, i.e., the signal reconstruction, corresponds to the inverse Fourier transform. In case of the Gabor basis, the reconstruction is sparse, i.e., only a finite number of function are necessary for reconstruction, contrary to the Fourier basis for which the reconstruction is continuous.

We evaluated the numerical cost of the OMP algorithm. For a dictionary size varying from 1000 × 50000 (N t = 1000; L = 20, K = 10, S = 250) to 1000 × 80000 (N t = 1000; L = 20, K = 10, S = 400), the computation time ranges approximately between 0.2 and 0.7 s on a standard personal computer (Intel(R) Xeon(R) CPU E5-2620 v2 at 2.10GHz, 32Gbytes memory, Dell, USA).

Orthogonal Matching Pursuit with complex dictionary

The application considered in this paper is the reconstruction the echoes from the periosteal and endosteal surfaces of the cortical bone. Accordingly, we restrain our analysis to the two first echoes, thus we set P = 2 in (1).

As a consequence, we expect only 2 non-zero components of x, which is much smaller than its dimension (M ). The determination of the two main echoes requires only two iterations using sparsity constrained OMP. Since the columns of the dictionary defined in (3) are complex-valued, the corresponding steps of the OMP have to be adapted. Each iteration step consists in searching the best matching column with the residual from the previous iteration. In the case of a complex dictionary, the matching between a particular complex column and the given signal r is [START_REF] Lu | Numerical implementation of matching pursuit for the analysis of complex ultrasonic signals[END_REF])

< r, d H >= d H r = Ce jφ (.5)
where r is the residual signal from last iteration, H denotes the complex conjugate transpose operator, C indicates the matching level and φ the phase of the residual.

The best correlated column from the dictionary D corresponds to the largest value of C. The corresponding position index in vector x is denoted

i max = max i {|D(:, i) H r|}, ( .6)
and the phase is given by φ = angle{D(:, i max ) H r}.

(.7)

It should be noticed that the residual cannot be iterated directly, since it is defined as the linear combination of two real Gaussian pulses instead of complex-values Gabor functions as introduced in (.1). After computation of (.6) and (.7), The real Gaussian pulses D r (i max ) are obtained using the relationship D r (:, i max ) = {D(:, i max )e -jφ }.

(.8)

The OMP processing is illustrated in Algorithm 1. Table 1: Ultrasound measurements: ∆t is the arrival time difference between the two echoes, ∆f is the shift of central frequency, σ f is related to the frequency bandwidth. The quality factor Q -1 11 is calculated from these three quantities and Ct.nBUA is calculated using the ultrasound measurements and Ct.

  Figure 1(a), while the two associated spectra S 1 (f ) and S 2 (f ) are shown in Figure 1(b). Those spectra were obtained by the temporal Fourier trans-

Figure 1 :

 1 Figure 1: An example of received echoes (for subject number 3) and illustration of the linear frequency dependence of the attenuation coefficient : (a) complete temporal ultrasound signal and retained segments s 1 (t) and s 2 (t); (b) frequency spectrum of the two echoes and (c) spectrum ratio ln(|S 2 (f )/S 1 (f )|) compared with the slope ∆f /σ 2 f obtained with Kuc's formula . The parameters t 1,2 , f 1,2 and σ t obtained with OMP in this study are indicated on the figures.

  using the same signal as for figure1.It can be observed in Figure2(b) that the retained parameters (best model number m) are associated with the two maxima, for all m, of the scalar product < s(t)|d(t; σ tm , f m , t m ))| > between the dictionary and the signal. The maximum of this scalar product can be interpreted as the maximum quality of the signal fit, as the dictionary vectors d are normalized. Accordingly, for each subject, we retained the 10 (out of 30) measurements with the highest values in order to remove the poorest measurements corresponding e.g., to a poor alignment of the probe with the bone surface. The final values of nBUA and Q -1 11 , and their errors, were obtained using the mean and the standard deviation calculated on the 10 retained measurements.

Figure 2 :

 2 Figure 2: Principle of the evaluation of the parameters (σ tm , f m , t m ) of the Gabor functions matching the two echoes using OMP. Original signal y(t) compared with the two reconstructed echoes s 1 (t) (blue) and s 2 (t) (black) (a) associated with dictionary indices m corresponding to the two largest values of the scalar product between the signal and the dictionary (b). (c): zoom around the first maximum.

  algorithm. The threshold used to discriminate cortical bone was set to one third of the apparent cortical bone density value (Dcort). Cortical thickness was calculated for each slice as the mean distance between the periosteal and endosteal contours. Finally, the mean and standard deviation of all slice cortical thickness was retained for each subject. Volumetric bone mineral density of the cortical bone in the ROI (Ct.vBMD) was obtained directly from Dcort. The mean and standard deviation of all slice Dcort was retained for each subject.

Figure 3 :

 3 Figure 3: A representative cross-sectional HR-pQCT slice of the radius for each subject. The region of interest corresponding the ultrasound measurement site is highlighted

Figure 4 :

 4 Figure 4: Correlation between quality factor Q -1 11 obtained by pulse echo and Ct.vBMD obtained by HR-pQCT.

Figure 5 :

 5 Figure 5: Correlation between Ct.nBUA and Ct.vBMD.

(

   which can be compared to our result). The differences between these equations may in part be due to the Ct.vBMD range which was different in the two studies: [1.15-1.65] g.cm -3 in(Sasso et al., 2008) compared [1.05- 1.22] g.cm -3 in the present study.

Algorithm 1

 1 Orthogonal Matching Pursuit algorithm Input: y ∈ R N , D ∈ C N ×M , target sparsity=2 Output: sparse solution x ∈ R M 1: Set I = (), r = y, D r = D(:, I) 2: for target sparsity do 3: i max = max i |D(:, i) H r| 4: I = (I, i max ) 5: φ = angle{D(:, i max ) H r} 6: D r (:, i max ) = {D(:, i max )e -jφ } 7: x(I) = D r (:, I) † y † denotes pseudo inverse 8: r = y -D r (:, I)x(I) 9: end for Tables

  to attenuation in cortical bone can be measured in vivo at the radius, which provide an information on bone quality as they were found to be highly correlated to Ct.vBMD values. Indeed, Ct.vBMD is an established biomarker of bone health related to bone porosity and mineralization[START_REF] Engelke | Quantitative Computed Tomography-Current Status and New Developments[END_REF]. If these results are confirmed in studies with a larger number and diversity of subjects, measurement of attenuation may be considered useful for assessing bone health. This can be combined with the measurement of cortical thickness, porosity and bulk wave velocities in multimodal cortical bone QUS evaluation methods. Future studies should investigate to which extent Ct.nBUA and Q -111 measured with the method of this study reflect the true (intrinsic) ultrasonic attenuation in cortical bone. Advanced signal processing techniques such as neural network

  2)where ζ Θ is a normalization parameter that ensures d(Θ, t) 2 = 1 (l 2norm). To build a suitable elementary atom of D, the columns use a discretization of the different parameters that characterize each Gabor function. The possible values of the set of parameters Θ are sampled on M discrete points. For this purpose, we define an a priori range for s, τ , and f , these ranges being subdivided into L, S and K regular intervals, respectively. Furthermore, each Gabor function is sampled atN t discrete time points t = [t 1 , t 2 , • • • , t Nt ].The dimension of the over-complete dictionary is

Table 2 :

 2 Th measured with Ht-pQCT. Cortical thickness (Ct.Th) and volumetric bone mineral density (Ct.vBMD) obtained from HR-pQCT images.

	Subject	∆t	∆f	σ f	Q -1 11	Ct.nBUA
	No.	(µs)	(MHz)	(MHz)	(%)	(dB.MHz -1 .cm -1 )
	1	1.67±0.02 0.21±0.15 0.94±0.03 5.1±3.4	3.7±2.5
	2	1.75±0.04 0.34±0.13 1.05±0.19 8.4±2.6	6.5±2.2
	3	1.70±0.02 0.46±0.07 0.81±0.02 13.5±3.1	9.5±2.3
	4	2.14±0.08 0.51±0.11 0.93±0.06 9.1±1.8	7.8±1.7
	5	2.32±0.04 0.37±0.16 0.89±0.06 6.8±3.6	5.8±3.1
	mean ± std 1.92±0.30 0.38±0.12 0.92±0.09 8.6±3.1	6.7±2.2
			HR-pQCT
	Subject No. Ct.Th (mm) Ct.vBMD (g.cm -3 )
	1	3.11±0.08	1.23±0.02
	2	3.13±0.02	1.13±0.02
	3	3.28±0.06	1.06±0.03
	4	3.43±0.02	1.14±0.02
	5	3.79±0.03	1.16±0.03
	mean ± std	3.35±0.28	1.14±0.06

Table 3 :

 3 Longitudinal bulk wave velocities calculated as v 11 = 2Ct.Th/∆t using values given in Tables1 and 2

	Subject No. v 11 (mm.µs -1 )
	1	3.72±0.14
	2	3.58±0.10
	3	3.86±0.12
	4	3.21±0.14
	5	3.27±0.08
	mean ± std	3.53±0.28

Table 4 :

 4 Normalized broadband ultrasonic attenuation (nBUA) values published in the scientific literature and in this study. Mean values are given, and standard deviation when available. * In Talmant et al. (2019), the reported attenuation normal to the direction of osteons is, at 4 MHz, 3.9 dB.cm -1 per percentage of porosity ; the value given in the table is calculated for a very moderate porosity of 5% characteristic of the bone of a healthy young adult.

					Frequency
	Reference		Species and skeletal sites	Number of samples	range or cen-ter frequency	nBUA MHz -1 cm -1 ) (dB
					(MHz)
	Lakes (1986)(Lakes et al., et al. 1986)	bovine femur (ex vivo)	1	1-7	∼3
	Lees and Klopholtz (1992)(Lees and Klopholz, 1992)	bovine femur (ex vivo)	4	0-30	∼4
	Han (1996)(Han et al., et al. 1996)	bovine femur (ex vivo)	5	0.3-0.7	5 -12
	Zheng al.(2007)(Zheng et al., 2007)	et	bovine femur (ex vivo)	8	2.25	4.91±0.65
	Sasso al.(2008)(Sasso et al., 2008)	et	bovine femur (ex vivo)	40	3.5-4.5	4.2±2.4
	Talmant et al. (2019) *	human femur (ex vivo)	35	2-8	4.9
			one-third dis-
	This paper		tal radius (in	5	3.5	6.7± 2.2
			vivo)

Acknowledgements

This work was supported by the GEP scholarship granted by the Guangzhou scholar and ANID ECOS200061 exchange programs. Jean-Gabriel Minonzio is supported by Grant ANID / FONDECYT / REGULAR / 1201311. The authors would like to thank Sylvie Fernandez and Christine Chappard (University Denis Diderot, CNRS, Osteo-Articular Bioengineering and Bioimaging (B2OA) 10, Avenue de Verdun, 75010 Paris, France) for the HR-pQCT measurements.

APPENDIX: Orthogonal Matching Pursuit

Dictionary

To implement the sparse signal processing method, the first issue is to discretize the Gabor functions into an over-complete dictionary D, where the columns of D are built resorting to the Gabor functions basis


Jan Behrends

Sthitadhi Roy

Michael H Kolodrubetz

Jens H Bardarson

Adolfo G Grushin

Landau levels, Bardeen polynomials and Fermi arcs in Weyl semimetals: the who's who of the chiral anomaly

Quantum anomalies describe the breaking of a classical symmetry by quantum fluctuations [START_REF] Bertlmann | Anomalies in quantum field theory[END_REF]. The chiral anomaly, the nonconservation of the chiral charge of three-dimensional Weyl fermions, is relevant to different domains in physics since Weyl fermions mediate the pion-decay into photons [START_REF] Stephen | Axial-Vector Vertex in Spinor Electrodynamics[END_REF][START_REF] Bell | A PCAC puzzle: π 0 → γγ in the σ-model[END_REF] and are emergent quasiparticles in Weyl semimetals [START_REF] Turner | Beyond band insulators: Topology of semimetals and interacting phases[END_REF][START_REF] Hosur | Recent developments in transport phenomena in Weyl semimetals[END_REF][START_REF] Miransky | Quantum field theory in a magnetic field: From quantum chromodynamics to graphene and Dirac semimetals[END_REF][START_REF] Jia | Weyl semimetals, Fermi arcs and chiral anomalies[END_REF][START_REF] Burkov | Topological semimetals[END_REF][START_REF] Armitage | Weyl and Dirac semimetals in three-dimensional solids[END_REF]. The physics is particularly transparent in the Landau level picture pioneered by Nielsen and Ninomiya [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF], requiring only basic quantum mechanics. In a magnetic field B, the conical Weyl dispersion evolves into Landau levels with a degeneracy proportional to |B| [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF]. Since momentum along the direction of the magnetic field remains a good quantum number, the Landau levels disperse in that direction, with the zeroth Landau level having a linear dispersion with a sign determined by the chirality; all other Landau levels have a quadratic dispersion. The zeroth Landau level of the left-and right-handed Weyl fermions furthermore connect at high energy. Consequently, an electric field (E) with a component along the dispersion transfers lefthanded fermions to right-handed fermions (or vice versa) resulting in a nonconservation of left and right particle numbers proportional to E • B [START_REF] Bertlmann | Anomalies in quantum field theory[END_REF][START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF][START_REF] Barry | Anomalies for pedestrians[END_REF].

Other fields, such as chiral pseudo-electromagnetic fields, torsion or curvature activate the chiral anomaly beyond E and B [START_REF] Parrikar | Torsion, parity-odd response, and anomalies in topological states[END_REF][START_REF] Sun | Chiral viscoelastic response in Weyl semimetals[END_REF][START_REF] Zubkov | Emergent gravity and chiral anomaly in dirac semimetals in the presence of dislocations[END_REF][START_REF] Cortijo | Elastic gauge fields in Weyl semimetals[END_REF][START_REF] Pikulin | Chiral Anomaly from Strain-Induced Gauge Fields in Dirac and Weyl Semimetals[END_REF][START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] You | Response properties of axion insulators and Weyl semimetals driven by screw dislocations and dynamical axion strings[END_REF][START_REF] Landsteiner | Gravitational Anomaly and Transport Phenomena[END_REF][START_REF] Gooth | Experimental signatures of the mixed axial-gravitational anomaly in the Weyl semimetal NbP[END_REF][START_REF] Ferreiros | Mixed axial-torsional anomaly in weyl semimetals[END_REF][START_REF] Lepori | Axial anomaly in multi-Weyl and triple-point semimetals[END_REF]. Weyl semimetals are ideal to probe the chiral anomaly in the presence of chiral pseudoelectromagnetic fields. To motivate this, recall that their low-energy degrees of freedom are pairs of chiral Weyl quasiparticles at topologically protected band touchings (Weyl nodes), separated in energy-momentum space by a four-vector b µ [START_REF] Armitage | Weyl and Dirac semimetals in three-dimensional solids[END_REF]. A space-and time-dependent b µ , as in strained or inhomogenously magnetized Weyl semimetals [START_REF] Cortijo | Elastic gauge fields in Weyl semimetals[END_REF] or Helium-3 [START_REF] Bevan | Momentum creation by vortices in superfluid 3 He as a model of primordial baryogenesis[END_REF][START_REF] Klinkhamer | Emergent CPT violation from the splitting of Fermi points[END_REF], generates chiral pseudomagnetic (B 5 = ∇×b) and pseudoelectric (E 5 = -∂ t b-∇b 0 ) fields, which couple with opposite signs to opposite chiralities [START_REF] Pikulin | Chiral Anomaly from Strain-Induced Gauge Fields in Dirac and Weyl Semimetals[END_REF][START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] Liu | Chiral gauge field and axial anomaly in a Weyl semimetal[END_REF][START_REF] Chernodub | Condensed matter realization of the axial magnetic effect[END_REF][START_REF] Schuster | Dissipationless conductance in a topological coaxial cable[END_REF][START_REF] Huang | Topological responses from chiral anomaly in multi-Weyl semimetals[END_REF]. These pseudofields enhance or generalize phenomena ranging from transport to interface physics [START_REF] Pikulin | Chiral Anomaly from Strain-Induced Gauge Fields in Dirac and Weyl Semimetals[END_REF][START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] Zhou | Topological invariants of metals and the related physical effects[END_REF][START_REF] Sumiyoshi | Torsional Chiral Magnetic Effect in a Weyl Semimetal with a Topological Defect[END_REF][START_REF] Cortijo | Strain-induced chiral magnetic effect in Weyl semimetals[END_REF][START_REF] Westström | Designer Curved-Space Geometry for Relativistic Fermions in Weyl Metamaterials[END_REF][START_REF] Gorbar | Pseudomagnetic lens as a valley and chirality splitter in Dirac and Weyl materials[END_REF][START_REF] Gorbar | Second-order chiral kinetic theory: Chiral magnetic and pseudomagnetic waves[END_REF][START_REF] Gorbar | Pseudomagnetic helicons[END_REF][START_REF] Roy | Tunable axial gauge fields in engineered Weyl semimetals: semiclassical analysis and optical lattice implementations[END_REF][START_REF] Shengyuan | Chirality-Dependent Hall Effect in Weyl Semimetals[END_REF][START_REF] Jiang | Topological imbert-fedorov shift in weyl semimetals[END_REF][START_REF] Cortijo | Emergent gravity in the cubic tight-binding model of Weyl semimetal in the presence of elastic deformations[END_REF]. Unlike B, which generates Landau levels dispersing in opposite directions for opposite chiralities, B 5 generates pseudo-Landau levels that disperse in the same direction for both chiralities [START_REF] Pikulin | Chiral Anomaly from Strain-Induced Gauge Fields in Dirac and Weyl Semimetals[END_REF][START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] Liu | Chiral gauge field and axial anomaly in a Weyl semimetal[END_REF][START_REF] Chernodub | Condensed matter realization of the axial magnetic effect[END_REF][START_REF] Landsteiner | Notes on anomaly induced transport[END_REF].

Applying the Landau level picture to chiral fields leads to puzzling conclusions: for example, because of the chirality independent dispersion of the zeroth pseudo-Landau level due to B 5 , E increases (or decreases) the number of fermions for both chiralities at a rate proportional to E • B 5 , giving an apparent nonconservation of the total charge. This is expressed as the so-called covariant anomaly [START_REF] Bertlmann | Anomalies in quantum field theory[END_REF][START_REF] Landsteiner | Notes on anomaly induced transport[END_REF]

∂ µ J µ 5,cov = 1 2π 2 (E • B + E 5 • B 5 ) , (1)
∂ µ J µ cov = 1 2π 2 (E 5 • B + B 5 • E) . (2)
Neither the covariant chiral (J µ 5,cov = J µ L,cov -J µ R,cov ) nor vector (J µ cov = J µ L,cov + J µ R,cov ) currents are conserved. In field theory, to explicitly restore charge conservation, the covariant currents are supplemented by Bardeen poly-nomials δj µ , which act like boundary conditions for the accumulated charge at the cut-off energy [START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF][START_REF] Gorbar | Consistent Chiral Kinetic Theory in Weyl Materials: Chiral Magnetic Plasmons[END_REF][START_REF] Bardeen | Anomalous Ward Identities in Spinor Field Theories[END_REF][START_REF] Gorbar | Origin of Bardeen-Zumino current in lattice models of Weyl semimetals[END_REF]. This procedure defines the consistent anomaly [START_REF] Bertlmann | Anomalies in quantum field theory[END_REF][START_REF] Landsteiner | Notes on anomaly induced transport[END_REF]

∂ µ J µ 5 = 1 2π 2 E • B + 1 3 E 5 • B 5 , (3)
∂ µ J µ =0. ( 4
)
The consistent anomaly conserves charge, and thus determines observables [START_REF] Gorbar | Consistent Chiral Kinetic Theory in Weyl Materials: Chiral Magnetic Plasmons[END_REF][START_REF] Gorbar | Chiral magnetic plasmons in anomalous relativistic matter[END_REF] by discarding unphysical responses [START_REF] Landsteiner | Notes on anomaly induced transport[END_REF][START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF][START_REF] Gorbar | Consistent Chiral Kinetic Theory in Weyl Materials: Chiral Magnetic Plasmons[END_REF][START_REF] Vazifeh | Electromagnetic Response of Weyl Semimetals[END_REF]. The Bardeen polynomials, however, sacrifice intuition of the covariant picture based on Landau levels and obfuscate the restoration of charge conservation in specific lattice implementations. Despite the field theory of the consistent and covariant anomalies being well understood for a long time [START_REF] Landsteiner | Notes on anomaly induced transport[END_REF], a simple physical picture of their origin on a lattice, with guaranteed charge conservation, is still missing. In this work we provide such a picture using as building blocks the Landau and pseudo-Landau levels. It leads to our two main results: first, we identify the Fermi arcs as a source for the covariant anomaly terms of Eq. ( 2) and relate them to the Bardeen polynomials. Second, we show that when B 5 > B the Fermi surface twists into a bowtie shape, a property central to our understanding of how the term E • B 5 redistributes charge within the sample. Our picture allows us to address as well the 1/3 disparity between the second term in Eq. ( 3) and Eq. ( 1). We argue that a necessary condition for its identification is a B 5 profile that spatially separates chiral charge creation and annihilation. Similar to how Landau levels simplified our understanding of the chiral anomaly [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF], we use pseudo-Landau levels (developed in Ref. 17) to provide a unified and simple lattice picture of the consistent and covariant anomalies with specific implications for experiment.

Our starting point is the Weyl semimetal model [START_REF] Vazifeh | Electromagnetic Response of Weyl Semimetals[END_REF] H

=v [sin(k y )σ x -sin(k x )σ y ] τ z + v sin(k z )τ y + mτ x + t i [1 -cos(k i )] τ x + v µ u µ b µ , (5)
with a = 1 the lattice constant, σ i (τ i ) spin (orbital) Pauli matrices, and ). Our results also apply to generalizations of Eq. ( 5) that model Dirac (e.g., Cd 3 As 2 , Na 3 Bi) and Weyl (e.g., TaAs family) materials [START_REF] Wang | Dirac semimetal and topological phase transitions in A3Bi (A =Na, K, Rb)[END_REF][START_REF] Wang | Three-dimensional Dirac semimetal and quantum transport in Cd3As2[END_REF][START_REF] Weng | Weyl Semimetal Phase in Noncentrosymmetric Transition-Metal Monophosphides[END_REF]. We further define γ µ = (τ x , iσ y τ y , -iσ x τ y , iτ z ) and γ 5 = iγ 0 γ 1 γ 2 γ 3 = σ z τ y , such that u µ = γ 0 γ µ γ 5 and the spacedependent chiral charge density is The charge density, ρ(x), is obtained by replacing γ 5 with the identity in Eq. [START_REF] Miransky | Quantum field theory in a magnetic field: From quantum chromodynamics to graphene and Dirac semimetals[END_REF]. While our lattice model [START_REF] Hosur | Recent developments in transport phenomena in Weyl semimetals[END_REF] includes a τ x term that explicitly breaks conservation of ρ 5 , we show in the Supplemental Material how this effect is controlled. The Hamiltonian (5) (with v = 1) derives from an effective field theory with action [START_REF] Jia | Weyl semimetals, Fermi arcs and chiral anomalies[END_REF] where ψ = ψ † γ 0 and repeated indices are summed. S yields two species of Weyl fermions of opposite chiralities, coupled to an external chiral field b µ and a vector field A µ , and separated by b µ for m = 0. The lattice regularization is given by the Wilson map k i → sin k i [START_REF]We restrict our analysis to isotropic Weyl nodes[END_REF], and m → m + t i (1 -cos k i ) [START_REF] Heinz | Lattice Gauge Theories: An Introduction[END_REF].

u µ = (σ z τ y , -σ x τ x , -σ y τ x , σ z ). For small b µ = (b 0 , b) and m 2 < v 2 |b 2 -b 2 0 |,
J 0 5 (x) ≡ ρ 5 (x) = n∈occ. ψ n (x)|γ 5 |ψ n (x) . ( 6
S = d 4 x ψ γ µ i∂ µ -A µ -b µ γ 5 -m ψ,
Spatial and temporal variations of b generate the chiral fields B 5 = ∇ × b and E 5 = -∂ 0 b. The simplest realization of B 5 occurs at the boundary of any Weyl semimetal with vacuum, where the Weyl node separation b µ goes to zero [START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] Chernodub | Condensed matter realization of the axial magnetic effect[END_REF][START_REF] Landsteiner | Notes on anomaly induced transport[END_REF]. For example, for a slab along y, b(y

) = b z [Θ(y -L/2) -Θ(y + L/2)]ẑ gives B 5 (y) = b z [δ(y -L/2) -δ(y + L/2)]
x, localized at the surface, which generates surface pseudo-Landau levels dispersing along ±k x , with opposite signs at each surface [START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF]. Their Fermi surface traces an arc, establishing the correspondence between surface pseudo-Landau levels induced by B 5 and topological surface states. Analogously, a uniform external magnetic field B = B ẑ parallel to the Weyl node separation leads to a spectrum hosting bulk Landau levels dispersing along ±k z , where the sign is set by the Weyl node chirality. When both B and a surface B 5 are present, Landau and pseudo-Landau levels coexist and the Fermi surface at the Fermi energy ε F (set to ε F = 0 hereafter) traces a rectangle [START_REF]The Fermi surface rectangle becomes a parallelogram for general orientations of B[END_REF], see Fig. 1(a) [START_REF] Bulmash | Quantum oscillations in Weyl and Dirac semimetal ultrathin films[END_REF][START_REF] Ominato | Magnetotransport in Weyl semimetals in the quantum limit: Role of topological surface states[END_REF].

The coexistence of Landau and pseudo-Landau levels provides an ideal platform to discuss the anomalies. Applying E = E z ẑ pumps charges of one chirality to the other, connecting Landau levels of B through the band bottom and realizing the E • B term in Eq. (3) [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF] [see Fig. 1(b)]. Similarly, since the pseudo-Landau levels disperse along ±k x on each surface, applying E = E x x depletes charges from one surface and generate charges on the other [Fig. 1(c)]. We can interpret [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF] this as an anomaly of each surface state due to E • B 5 . In contrast to the chiral anomaly in the absence of chiral fields, where the total charge is locally conserved, the spatial separation of the two surfaces leads to an apparent violation of local charge conservation, as in Eq. ( 2). Our pseudo-Landau level picture demonstrates the surface origin of the covariant anomaly, ∼ E • B 5 .

Our picture tracks how charge is explicitly conserved. The spectral flow between the pseudo-Landau levels at each surface happens via the bulk Landau levels connecting them, fixing local charge conservation. The effect of E results in an adiabatic shift of k x , which via positionmomentum locking [START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF] generates a Hall current

δj = - 1 2π 2 b z E x ŷ. ( 8
)
We can interpret this as the net current flowing along ŷ between surfaces through the bottom of the band. In the bulk, b z is constant leading to ∇ • δj = 0 and no accumulation of charge. At the surface, the Weyl node separation varies, leading to a finite divergence of the spatial current, positive on one surface and negative on the other. Similarly, more general profiles of B and B 5 can be understood in terms of Landau and pseudo-Landau levels. For instance, a uniform bulk By construction, the corresponding current J µ cov misses information from states away from the Fermi level, and thus it is not conserved as dictated by Eq. ( 2), cf. Supplemental Material. The consistent current, J µ , is obtained from the covariant current by adding the Bardeen polynomials δj µ such that ∂ µ (J µ cov +δj µ ) = ∂ µ J µ = 0 [START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF][START_REF] Gorbar | Consistent Chiral Kinetic Theory in Weyl Materials: Chiral Magnetic Plasmons[END_REF]. Using Eq. ( 2) and the definition of the pseudo-fields,

δj 0 = 1 2π 2 b • B, δj = 1 2π 2 (b 0 B -b × E) . (9)
Comparing ( 9) to (8) of our first example (Fig. 1), we identify the latter as a part of the Bardeen polynomials [START_REF] Landsteiner | Notes on anomaly induced transport[END_REF]. The benefit of the Landau level approach is its intuitive interpretation: in the first example, E 5 = 0 and the finite E • B 5 pumps charge from one surface to another via the anomalous Hall effect [Eq. ( 8)] through

∂ t J 0 cov = E • B 5 /(2π 2 ).
Our second example, Fig. 2, can be interpreted similarly. In the bulk charge grows as Moving on to spatio-temporally varying Weyl node separations, consider first b = -E 5 tẑ, yielding uniform E 5 ẑ. A B = 0 activates the second term of (2), suggesting that charge is created at the Fermi surface at a rate E 5 • B/(2π 2 ). For our lattice model ( 5), E 5 shifts the band bottom, pushing charge above a fixed energy [Fig. 3(a)]. Rigid shifts of the band conserve total charge, giving the consistent picture of the anomaly. However, if one insists on only considering the low-energy gray region in Fig. 3(a), the charge appears to emerge from the vacuum, as expected from the covariant anomaly. The connection between these pictures is shown in Fig. 3(b), where charge growth near the Fermi surface equals charge loss near the band bottom, which in turn equals the Bardeen polynomial δj 0 , Eq. ( 9).

∂ t J 0 cov = E • B 5 /( 2π
We end by addressing the factor of 1/3 disparity between the prefactors of the E 5 • B 5 anomalies. This difference implies that the band bottom current must add -2/3 to the Fermi surface contribution, irrespective of the precise pseudo-field profile. One may argue that this factor arises from the topological nature of the Bardeen polynomials [START_REF] Gorbar | Origin of Bardeen-Zumino current in lattice models of Weyl semimetals[END_REF], yet we find that topology alone does not explain the conditions which give 1/3 for a generic lattice model.

To illustrate the conditions for isolating the 1/3, consider [START_REF] Roy | Tunable axial gauge fields in engineered Weyl semimetals: semiclassical analysis and optical lattice implementations[END_REF] ]. This profile maintains the useful property that the two chiralities (eigenstates of γ 5 ) are well separated in momentum space [Fig. 4(c)]. Adding E 5 B 5 to produce a chiral anomaly, chiral charge ( 6) is created/annihilated in spatially separated regions [Fig. 4(d)], at a rate that closely follows the spatial profile of B 5 as expected from Eq. (3) [Fig. 4(e)]. The Fermi surface contribution to the chiral charge is calculated as

b = c B 2 erf y -L 4 √ 2ξ -erf y + L 4 √ 2ξ -1 x + b z ẑ, (10)
dρ F S 5 (y) db = n ψ n |γ 5 Π y |ψ n ψ n |∂ b H|ψ n δ(µ -ε n ).
(11) with Π y the projector to the position y, cf. Supplemental Material. Comparing the chiral charge pumped at the Fermi surface and that in the full spectrum to the expected value proportional to B 5 , we see that they are close, but neither perfectly reproduces the field theory prediction.

To determine the factors behind this mismatch we note that the erf profile minimizes the spatial overlap between regions of chiral charge creation and annihilation (∂ρ 5 /∂b z > 0 and < 0 at y ≈ +L/4 and y ≈ -L/4 respectively). However, for both chiralities, there exist regions in momentum space [light/dark regions in Fig. 4(b)] where the states are not well-localized in real space. These regions spread out with increasing ξ, such that the erf profile evolves into a sinusoid-like function, resulting in a poorer match between the B 5 (y) profile and the chiral charge response, cf. Supplemental Material. The sensitivity of the result to these regions, highlighted by the failure of the sinusoid-like profile where no clear spatial separation exists, implies that the value 1/3 is modified by such effects for generic profiles of B 5 . This detrimental overlap is minimized if L ξ a is satisfied, allowing the intriguing possibility that an exact 1/3 may be recovered in this limit. Finally, all anomaly terms present finite size and quadratic corrections to the low-energy field theory [START_REF] Jia | Weyl semimetals, Fermi arcs and chiral anomalies[END_REF], which we discuss in the Supplemental Material.

In summary, we have provided an intuitive lattice picture that is based on Landau and pseudo-Landau levels and connects the covariant and consistent anomalies. We explicitly identified the Hall current as the Bardeen polynomial that connects the covariant anomalies of the Fermi arcs and restores charge conservation, most notably when the Fermi surface knots into a bowtie. We expect that the bowtie Fermi surface and its response to external fields will be important to understand strained Weyl semimetals experimentally, in particular their transport properties.

Our work highlights that measuring the consistent or covariant anomaly (e.g., the factor of 1/3) depends on whether the experimental probe is sensible to only the Fermi surface or rather the entire Fermi sea. Additionally, perturbations such as strain allow other model parameters to depend on position, e.g, the Fermi velocity, as well as additional terms in Eq. ( 7) [START_REF] Arjona | Rotational strain in Weyl semimetals: A continuum approach[END_REF]. Our work motivates the study of these questions to interpret incipient experiments in strained Weyl semimetals.

For completeness, we provide a summary of the relevant facts of consistent and covariant anomalies derived from quantum field theory. As discussed in the main text electric and magnetic fields that satisfy E•B = 0 generate an anomalous imbalance between left and right movers. Additionally pseudo-magnetic (B 5 ) and pseudo-electric (E 5 ) field-the fields that couple with opposite signs to opposite chiralities-will result in the same chiral imbalance, since the contributions of the two chiralities to the chiral current J µ 5,cov = J µ L,cov -J µ R,cov (µ = 0, 1, 2, 3) will add up. However, these considerations imply also that the vector current

J µ cov = J µ L,cov +J µ R,cov is not conserved. Mathematically ∂ µ J µ 5,cov = 1 2π 2 (E • B + E 5 • B 5 ) , (12a)
∂ µ J µ cov = 1 2π 2 (E • B 5 + E 5 • B) , (12b)
a result which can be obtained either diagrammatically [START_REF] Bertlmann | Anomalies in quantum field theory[END_REF][START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF] or semiclassically [START_REF] Pikulin | Chiral Anomaly from Strain-Induced Gauge Fields in Dirac and Weyl Semimetals[END_REF][START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] Liu | Chiral gauge field and axial anomaly in a Weyl semimetal[END_REF][START_REF] Gorbar | Consistent Chiral Kinetic Theory in Weyl Materials: Chiral Magnetic Plasmons[END_REF]. Eqs. ( 12) define the covariant chiral anomaly [START_REF] Landsteiner | Notes on anomaly induced transport[END_REF][START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF] of the covariant vector and axial currents J µ cov and J µ 5,cov respectively. This result is troublesome: a theory which does not conserve charge is unphysical, yet the fields B 5 and E 5 can be physically realized in the solid state [START_REF] Pikulin | Chiral Anomaly from Strain-Induced Gauge Fields in Dirac and Weyl Semimetals[END_REF][START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF][START_REF] Liu | Chiral gauge field and axial anomaly in a Weyl semimetal[END_REF]. At the quantum field theory level curing the nonconservation of charge amounts to acknowledging that Eqs. [START_REF] Parrikar | Torsion, parity-odd response, and anomalies in topological states[END_REF] are determined solely by the Fermi surface, but that counteracting currents are necessary at energies comparable to the cut-off energy [START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF]. To ensure this constraint is satisfied, two additional Chern-Simons currents, δj µ and δj µ 5 , known as Bardeen polynomials, are added, which transfer the necessary charge and guarantee charge conservation [START_REF] Bardeen | Anomalous Ward Identities in Spinor Field Theories[END_REF]. They relate the covariant currents (J µ cov , J µ 5,cov ) to the physical consistent currents (J µ , J µ 5 ) through J µ 5 = J µ 5,cov + δj µ 5 and J µ = J µ cov + δj µ , respectively, and take the form

δj 0 = 1 2π 2 A 5 • B; δj = 1 2π 2 A 0 5 B -A 5 × E ; ( 13
)
δj 0 5 = 1 6π 2 A 5 • B 5 ; δj 5 = 1 6π 2 A 0 5 B 5 -A 5 × E 5 , (14)
where the chiral gauge field A µ 5 defines the pseudo-fields B 5 = ∇ × A 5 and E 5 = -∇A 0 5 -∂ 0 A 5 . Together with Eqs. [START_REF] Parrikar | Torsion, parity-odd response, and anomalies in topological states[END_REF] results in the consistent anomaly [START_REF] Landsteiner | Anomalous transport of Weyl fermions in Weyl semimetals[END_REF][START_REF] Gorbar | Consistent Chiral Kinetic Theory in Weyl Materials: Chiral Magnetic Plasmons[END_REF]

∂ µ J µ 5 = 1 2π 2 E • B + 1 3 E 5 • B 5 ; ∂ µ J µ = 0. ( 15
)
Eqs. ( 12) and ( 15) summarize the two versions of the quantum anomalies.

We note that the consistent and covariant currents have been discussed in the context of the two-dimensional quantum Hall effect [START_REF] Stone | Gravitational anomalies and thermal Hall effect in topological insulators[END_REF]. In this case the edge current is the sum of the edge consistent current plus the anomaly inflow from the 2+1D bulk Chern Simons term, which combined form the covariant edge current [START_REF] Chandrasekharan | Anomaly cancellation in 2+1 dimensions in the presence of a domain wall mass[END_REF]. In a similar way, consistent and covariant currents can be defined for a 3+1D Weyl fermion living at the edge of a 4+1D quantum Hall effect.

Details of the lattice calculations

We are interested in the response of a Weyl semimetal to external fields E and E 5 . A time-dependent vector potential A = -Et gives an electric field, and analogously a time-dependent node separation b = -E 5 t. We express the time-derivative of any expectation value f in terms of the response to changes in A and b,

df dt = -E • ∂ A f -E 5 • ∂ b f, (16)
which is especially handy when evaluating timederivatives proportional to E and E 5 , such as ∂ µ j µ and ∂ µ j µ 5 where f corresponds to the densities ρ or ρ 5 .

The response of the density to external fields E and E 5 is the response to a shift A → A + dA and b → b + db respectively. Explicitly, the chiral density at site y reads

ρ 5 (y) = n∈occ. ψ n |γ 5 Π y |ψ n , (17)
where Π y projects onto site y. When performing an infinitesimal shift of the Weyl node separation b → b + db, the single-particle states change as

|ψ n (b + db) = |ψ n + db • m =n |ψ m ψ n |∂ b H|ψ m ε n -ε m (18)
resulting in a change of the chiral density

dρ 5 (y) db = n∈occ. m =n ψ n |γ 5 Π y |ψ m ψ n |∂ b H|ψ m ε n -ε m + h.c. ( 19
) Note that the sum goes over the states n that are initially occupied at a node separation b.

The Fermi surface contribution, on the other hand, just counts the subset of states that are lifted above the Fermi level when shifting b → b + db, minus the states that are pushed below the Fermi level, dρ F S 5 = dρ + 5 -dρ - 5 with

dρ + 5 (y) = εn<µ, εn+dεn>µ ψ n |γ 5 Π y |ψ n (20) dρ - 5 (y) = εn>µ, εn+dεn<µ ψ n |γ 5 Π y |ψ n (21)
or, expressed in terms of Heaviside functions,

dρ ± 5 (y) = ± n ψ n |γ 5 Π y |ψ n Θ(±dε n ) × [Θ(µ -ε n ) -Θ(µ -(ε n + dε n ))] . (22)
To first order in perturbation theory, the infinitesimal change in the energies is

dε n = db • ψ n |∂ b H|ψ n (23)
and the sum of Heaviside functions reduces to

Θ(µ -ε n ) -Θ(µ -(ε n + dε n )) = dε n δ(µ -ε n ), (24)
which finally gives

dρ F S 5 (y) db = n ψ n |γ 5 Π y |ψ n ψ n |∂ b H|ψ n δ(µ -ε n ) (25
) as used in the main text. Numerically, we use a Lorentz approximation to the delta function, δ(x) = lim η→0 η/(π(η 2 + x 2 )).

The quest for 1/3

Two effects challenge the calculation of the chiral anomaly due to E 5 • B 5 on a lattice. First, B 5 always averages to zero over the whole sample [START_REF] Grushin | Inhomogeneous Weyl and Dirac Semimetals: Transport in Axial Magnetic Fields and Fermi Arc Surface States from Pseudo-Landau Levels[END_REF], which implies that a length scale ξ exists that characterizes the typical size of the variations in B 5 . Second, a mass term m changes the response of the chiral charge ∂ µ j µ 5 .

Effect of ξ

We demonstrate the lattice effect set by the length ξ explicitly by using the same profile of the node separation b as in the main text [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF]. In Fig. 5 (a) and(b), we show how the responses of the total chiral charge ρ 5 and the charge around the Fermi surface ρ FS 5 are affected by different sizes of ξ respectively. The two conditions for ξ discussed in the main text that allow for the identification of 1/3 are both demonstrated explicitly. When ξ ≈ a the variation of B 5 occurs on the scale of the lattice constant forcing the density variations Eqs. [START_REF] Landsteiner | Gravitational Anomaly and Transport Phenomena[END_REF] and [START_REF] Liu | Chiral gauge field and axial anomaly in a Weyl semimetal[END_REF] to interpolate smoothly between lattice sites such that the quantum field theory predictions breaks down. When ξ L is not satisfied, tails of the regions with positive and negative B 5 (y) (positive and negative chiral charge creation) overlap, leading to notable finite size effect and again conflicting with the field theory predictions. 5)] with periodic boundary conditions in y-directions of length L = 360 and cB = 1.25, cf. Eq. ( 10). We show results for different widths ξ of the regions with nonzero B5. The best matching of the field theory predictions are found for ξ = 6, which mostly satisfies the condition a ξ L. For ξ = 2, the condition a ξ is not satisfied, making lattice scale effects extremely relevant and deviating from the field theory limit. For ξ = 30, tails of the two chiral charge pumping regions centered at ±L/4 overlap, leading to notable finite size effects.

Effect of a finite mass

The response of the chiral and total charge to fields E and E 5 computed on a lattice can deviate from the expectation based on the simplest quantum field theory due to corrections set by the mass term, as mentioned in the main text. In this section, we carefully investigate one example, the electromagnetic contribution to the chiral anomaly on a lattice, ∂ µ j µ 5 = E • B/(2π 2 ) and argue how a mass term influences the chiral anomaly. We further show that the mass term does not play the same role for the covariant anomaly ∂ µ J µ cov = E 5 • B/(2π 2 ). The main source of lattice corrections to the anomaly is that the Wilson fermion Hamiltonian [START_REF] Kenneth | Confinement of quarks[END_REF] used for all tight-binding calculations, Eq. ( 5), has a momentumdependent mass term

M k = m + t i (1 -cos k i ), (26)
that ensures the absence of copies of the Weyl nodes, or doublers, on the lattice [START_REF] Kenneth | Confinement of quarks[END_REF]. When t > 0, the minimum of M k is at the Γ point and its maximum at (π, π, π); for m = 0, as set in the main text, the term is zero at k = 0 and increases away from Γ. The term M k sets the momentum-independent mass term m in the action (7) that couples both chiralities and changes the response to electric and magnetic fields [START_REF] Jackiw | Radiatively Induced Lorentz and CPT Violation in Electrodynamics[END_REF][START_REF] Pérez-Victoria | Exact Calculation of the Radiatively Induced Lorentz and CPT Violation in QED[END_REF][START_REF] Grushin | Consequences of a condensed matter realization of Lorentz-violating QED in Weyl semimetals[END_REF]. Such a term results in an additional classical contribution to the chiral anomaly, ∂ µ j µ 5,class = 2mi ψγ 5 ψ [START_REF] Bertlmann | Anomalies in quantum field theory[END_REF], that is zero in all equilibrium situations we consider, which we confirm numerically. Instead, another consequence of m affects the chiral anomaly: when m = 0 in the field theory, the eigenfunctions of the corresponding Hamiltonian are no longer eigenfunctions of the chiral matrix γ 5 . Instead, close to the Weyl node of chirality χ, the overlap with γ 5 for the low-energy eigenstates with a linear dispersion can be obtained analytically to be with b 2 = b µ b µ , i.e., the chiral density around the Weyl is reduced by a nonzero mass term. We note that a similar factor appears in the effective action for this model and is related to the Weyl node separation [START_REF] Pérez-Victoria | Exact Calculation of the Radiatively Induced Lorentz and CPT Violation in QED[END_REF][START_REF] Grushin | Consequences of a condensed matter realization of Lorentz-violating QED in Weyl semimetals[END_REF].

ψ χ |γ 5 |ψ χ = χ 1 - m 2 |b 2 | , (27) (a) (b) (c) (d) (e)
On a lattice, the eigenstates of the Hamiltonian (5) are only eigenstates of γ 5 for those momenta k where M k = 0. In particular, the overlap with γ 5 of states of chirality χ = ± close to the Weyl nodes at χk W is

ψ χ χk W |γ 5 |ψ χ χk W = χ 1 - M 2 k W |b 2 | , (28)
i.e., it scales with the value of M k at the Weyl nodes, m ≡ M k W . In the presence of m = 0 the eigenstates are not strictly eigenstates of γ 5 yet in order to define the chiral anomaly using Eq. ( 6) we effectively assume that they are. An alternative definition of chirality that does not possess this problem is to partition the Brillouin zone into two and declare left and right chiral charges as done in Ref. [START_REF] Roy | Chern numbers and chiral anomalies in Weyl butterflies[END_REF]. While it has the advantage of having a clear definition of chirality, it does not reduce to any representative matrix at the field theory level and so we retain the first definition. With these definitions it is tempting to upgrade the chiral anomaly to

∂ µ j µ 5 = 1 2π 2 1 - m2 b 2 E • B. (29)
It is interesting to note that a similar equation would be obtained if one interprets the prefactor as an effective renormalization of the chiral electric charge e 5 = e 1 -m2 /b 2 that would enter the chiral anomaly throughout the coupling to an external chiral gauge field e 5 A µ,5 in analogy with ordinary charge [START_REF] We | Juan for this observation[END_REF]. We have unsuccessfully attempted to derive Eq. ( 29) from the effective field theory responses for the action Eq. ( 7) discussed in [START_REF] Jackiw | Radiatively Induced Lorentz and CPT Violation in Electrodynamics[END_REF][START_REF] Pérez-Victoria | Exact Calculation of the Radiatively Induced Lorentz and CPT Violation in QED[END_REF][START_REF] Grushin | Consequences of a condensed matter realization of Lorentz-violating QED in Weyl semimetals[END_REF]. However, our numerical results on the corrections presented below suggest a richer structure when the mass is non-zero, which has motivated us to leave the precise connection between the quantum field theory discussed in these works and our results for a separate study. At the end of the section we specify the conditions where our numerics coincide with the massless field theory expectation validating our results in the main text, up to a given order in momentum.

To gain insight on the above expectation numerically, we investigate how the chiral anomaly changes with m. Differentiating between the total charge and the charge around the Fermi surface, we compute the response to external fields E and E 5 via Eqs. ( 19) and ( 25) on a lattice with periodic boundary conditions and a constant B.

We show the response of the chiral charge ρ 5 in Fig. 6(a) and (b). In panel (a), we show the convergence with system size, whereas in panel (b), we show that, in the limit of L → ∞, the response of the chiral charge scales with 1 -m2 /b 2 as dictated by Eq. ( 29).

Taking into account just states around the Fermi surface via Eq. ( 25), the chiral charge ρ FS 5 responds approximately in the same way as shown in Fig. 6(b), cf. panel (c). Thus, the consistent and covariant anomalies are the same for the term E • B, even in the presence of a mass m. In panel (d), we show the response of the total charge around the Fermi surface to a field E 5 : due to the covariant anomaly, ∂ρ FS /∂A 5 = 0. In contrast to ∂ρ FS 5 /∂A, the contribution ∂ρ FS /∂A 5 is not much affected by m = 0.

The chiral-symmetry breaking term M k has more intricate consequences for the anomaly in the chiral charge due to E 5 • B 5 . We find numerically that the dependence of ∂ρ 5 /∂b z and ∂ρ FS 5 /∂b z on m are different. In Fig. 7(a), we show the response of the total chiral charge for different values of m with the profile of b according to Eq. [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF]. It allows us to obtain panel (b), where the value of the integral over the curves in Fig. 7(a) as a function of m is compared to the result expected from the massless quantum field theory, i.e. the 1/3 prefactor difference between the second term of equations ( 13) and [START_REF] Cortijo | Elastic gauge fields in Weyl semimetals[END_REF] shown as a dashed line. ∂ρ 5 /∂b z is closest to the massless field theory prediction for a value m ∼ 0.1v, close to m = 0. This contrasts with the E • B analyzed above where the massless field theory expectation was reached when m = 0. Finally, the effect of m on ∂ρ FS 5 /∂b z is shown in Fig. 7(c) and (d) for various values of the level broadening η. The integral over one peak decreases with m, as visible in Fig. 7(e), different from the total charge contribution.

In conclusion, the massless quantum field theory expectation is recovered numerically when m ∼ m ∼ 0. This is satisfied when quadratic corrections are negligible and we choose m ∼ 0, a regime which is largely satisfied by our parameter choices in the main text. The small yet rich deviations that we have quantified numerically (Figs. 6 and7) indicate that further effort beyond the scope of this work is required to precisely connect with quantum field theory expectations with a finite mass.

3 .

 3 the model has one pair of Weyl nodes near Γ [46, 47]. Unless stated otherwise, we set m = 0 and t = 2v/ √ When b is oriented along a reciprocal lattice vector, this parameter choice gives two Weyl nodes located at ±b[1 + O(b 2 0 )] + O(b 5 j ) and energies ±vb 0 [1+O(b 2 j )]+O(b 3 0

)FIG. 1 .

 1 FIG. 1. Anomalies due to E•B and E•B5. (a) The spectrum of Hamiltonian (5) for a slab finite along y, showing the coexisting Fermi arcs and Landau levels of B = B ẑ. The color scale denotes the wave function's position along y. (b), (c) The Landau (pseudo-Landau) levels of B (B5) on the red (blue) triangle plane show the anomaly in the presence of E = Ez ẑ (E = Ex x), where the spectral flow at the Fermi surface is shown by green (orange) arrows. We implement E using the gauge choice k → k -Et. The left panels show the occupation at t = 0, while the right panels show the dispersion for t > 0 with Ex,zt = 0.08 and the same momenta occupied as in the left panels, cf. Supplemental Material. Here L = 100 and B = 1/ √ B = 11.2.

  B 5 ẑ arises from b = B 5 yx. Its spectrum [Fig. 2(a)] shows a characteristic butterfly Fermi surface, obtained from the rectangular Fermi surface Fig. 1(a) by noticing that the two bulk Landau levels have lengths B ± B 5 . When B 5 > B, the Fermi surface twists, leading to Fig. 2(a) (found to lead to peculiar quantum oscillation signals in an unrelated context[START_REF] Bovenzi | Twisted Fermi surface of a thin-film Weyl semimetal[END_REF]). An E = E ẑ parallel to B 5 [Figs.2(b) and (c)] makes the bulk gain charge above the Fermi level (upward arrows), while the surface loses charge (downward arrows). This is consistent with the lattice numerics [Fig.2(e)] where the spatial profile of the charge relative to that at t = 0 is shown for different times.Our previous examples (Figs.1 and 2) are a consequence of the covariant anomaly which only considers the depletion and growth of charges at the Fermi level.

FIG. 2 .

 2 FIG. 2. Anomaly due to E • B5 with a constant bulk B5 = B5 ẑ. (a) Energy spectrum with the same color coding and boundary conditions as Fig. 1. (b) The pseudo-Landau levels of B5 show an anomaly when E = E ẑ with the spectral flow indicated by the green arrows. The left panel shows the occupation at kx = 0.2 [gray plane in (a)] at t = 0; in the right panel Et = 0.05. (c) Left: The spatial profile of the charge δρ relative to that at t = 0, normalized by the total charge ρ0 for different times (colors), showing a surface to bulk charge redistribution. The y axis is chosen such that only the (positive) bulk contribution is visible, which is exactly canceled by the (negative) surface contribution. Right: Derivative of the density with respect to the vector potential A for different times, which equals ∂A z ρ = B5/(2π 2 ) in the bulk (dashed line). The derivatives with respect to vector potential and time are related via ∂tρ = -E • ∂ A ρ. In (a) and (b), L = 100 and 5 = 1/ √ B5 = 11.1; in (c), L = 200 and 5 = 15.8.

FIG. 4 .

 4 FIG. 4. Anomaly due to E5 • B5. (a) The Weyl node separation (10) and its corresponding B5. (b) Energy spectrum with a real space color encoding that reflects the periodic boundary conditions. (c), (d) Where B5 > 0 [red in (d)], a parallel E5 pushes left-and right-handed chiral charges above and below the Fermi level respectively, annihilating chiral charge. Where B5 < 0 (blue), chiral charge is created. (e) The spatial distribution of chiral charge creation and annihilation follows and equals B5 (red solid line) for charges ρ FS 5 traversing the original Fermi surface (blue dashed line). With all bands, the total chiral charge creation is ∼ 1/3 of that at the Fermi surface, as predicted by the consistent anomaly (black dotted line). We use L = 360, a B5 broadening ξ = 6 and peak height cB = 1.25, and a delta function broadening η = 10 -6 v, cf. Supplemental Material.

  2 ) locally while at the surface charge is depleted since B 5 has the opposite sign. The corresponding current[START_REF] Armitage | Weyl and Dirac semimetals in three-dimensional solids[END_REF] pumps charge from the surface to the bulk. Since b = B 5 yx, there are local currents in the bulk δj ∼ B 5 Ey ŷ and ∇ • δj = B 5 • E/(2π 2 ) is precisely the growth rate of local charge [Fig. 2(c)], reconciling the Fermi surface (covariant) picture with the charge conserving (consistent) picture on the lattice [57].

  with c B a constant and erf(x) the error function [Fig. 4(a)]. In a periodic system of length L in the y direction, this profile realizes interfaces at ±L/4 between two Weyl semimetals with different node separations, connected by regions where b smoothly changes to give B 5 = B 5 (y)ẑ with B 5 (y) = c B √ 2πξ [e -(y-L 4 ) 2 /(2ξ 2 )e -(y+ L 4 ) 2 /(2ξ 2 )

FIG. 5 .

 5 FIG. 5. (a) Responses ∂ρ5/∂bz and (b) ∂ρ FS 5 /∂bz for our lattice model [Eq. (5)] with periodic boundary conditions in y-directions of length L = 360 and cB = 1.25, cf. Eq. (10). We show results for different widths ξ of the regions with nonzero B5. The best matching of the field theory predictions are found for ξ = 6, which mostly satisfies the condition a ξ L. For ξ = 2, the condition a ξ is not satisfied, making lattice scale effects extremely relevant and deviating from the field theory limit. For ξ = 30, tails of the two chiral charge pumping regions centered at ±L/4 overlap, leading to notable finite size effects.

FIG. 6 .

 6 FIG. 6. (a) Convergence of the response ∂ρ5/∂Az for a tightbinding system governed by the Hamiltonian (5) with periodic boundary conditions in y-directions of length L and a node separation b = bẑ. A magnetic field B = B ẑ that satisfies periodic boundary conditions is included, with a magnetic length B = √ 2πL. The crosses denote the numerical results from tight-binding calculations , the solid line the results of a fit with f (L) = c0 + c1 exp(-L/ξ), and the dashed line the coefficient c0, i.e., the result in the limit L → ∞. The color denotes different values of m, as shown in panel (b), and we further choose b = 0.4. (b) The response ∂ρ5/∂Az in the limit L → ∞ as a function of m2 /b 2 . The different symbols denote different combinations of b and t and the gray line shows the expectation from the lattice model (29). (c) Response of the chiral charge around the Fermi surface to E ∂ρ FS 5 /∂A, as a function of the level broadening η for b = 0.4. The dashed lines denote the result from panel (a), limL→∞ ∂ρ5/∂A. (d) Response of the total charge around the Fermi surface to an axial field, ∂ρ FS /∂A5, as a function of the level broadening η. The dashed line denotes the field theory result ∂ρ/∂A5 = B/(2π 2 ) as a guide for the eyes.

FIG. 7 .

 7 FIG. 7. (a) Response ∂ρ5/∂bz for a tight-binding system governed by the Hamiltonian (5) with periodic boundary conditions in y-directions of length L = 360. An axial field B5 = B5(y)ẑ generated by b(y), Eq. (10), with ξ = 6 and cB = 1 that satisfies periodic boundary conditions is included. The different colors denote different values of m (cf. panel (b)) and the bold red line shows the profile of B5(y). (b) The integral over y around the peak is computed numerically for different values of m. The dashed line denotes the value of the integral that is expected from quantum field theory. In (c) and (d), we show the response of the chiral charge around the Fermi surface ∂ρ FS 5 /∂bz as a function of the level broadening η with the different colors denoting different values of η. The bold red line shows the profile of B5(y), the filled region is minimal difference between B5 and ∂ρ FS 5 /∂bz. In (e), we show the numerically computed integral over y around the peak in B5 for different values of m. The dashed line denotes the value of the integral that is expected from quantum field theory.

Supplemental MaterialSummary of consistent and covariant anomalies

Acknowledgements.-We thank B. Bradlyn, Y. Ferreiros, R. Ilan, F. de Juan, K. Landsteiner and D. Pesin for insightful discussions, and especially N.G. Goldman for discussions that initiated this work and comments on the manuscript. This research was funded by the Marie Curie Program under EC Grant Agreement No. 653846 (A.G.G.), the ERC Starting Grant No. 679722 and the Knut and Alice Wallenberg Foundation 2013-0093 (J.H.B.), the EPSRC Grant No. EP/N01930X/1 (S.R.), and supported by the Laboratory Directed Research and Development (LDRD) funding from Lawrence Berkeley National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DEAC02-05CH11231, and the DOE Basic Energy Sciences (BES) TIMES initiative as well as UTD Research Enhancement Funds (M.H.K.).