Multi-level Monte Carlo computation of the hadronic vacuum polarization contribution to ( g μ −2)
Mattia Dalla Brida, Leonardo Giusti, Tim Harris, Michele Pepe
MMulti-level Monte Carlo computation of the hadronic vacuumpolarization contribution to ( g µ − Mattia Dalla Brida a,b , Leonardo Giusti a,b , Tim Harris a,b , Michele Pepe b a Dipartimento di Fisica, Universit`a di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy b INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
The hadronic contribution to the muon anomalous magnetic moment a µ = ( g µ − / a µ .As the approach is applicable to other computations affected by a signal-to-noise ratio problem, ithas the potential to unlock many open problems for the nuclear and particle physics community. I. INTRODUCTION
The current experimental value of the muon anoma-lous magnetic moment a µ = 11659208 . . × − bythe E821 experiment has the remarkable precision of 0 . .
14 ppm by the end of its operation [2].The Standard Model (SM) prediction includes contribu-tions from five-loop Quantum Electrodynamics, two-loopWeak interactions, the Hadronic leading-order VacuumPolarization (HVP) and (the much smaller) HadronicLight-by-Light scattering (HLbL), see Ref. [3] and ref-erences therein. The overall theoretical uncertainty isdominated by the hadronic part. So far , lacking pre-cise purely theoretical computations, the hadronic con-tributions have been extracted (by assuming the SM)from experimental data via dispersive integrals (HVP& HLbL) and low-energy effective models supplementedwith the operator product expansion (HLbL). This leadsto a µ = 11659181 . . × − (0.37 ppm) [3], whichdeviates by 3 − .
6% to roughly 2% [5–11] corresponding to an overallerror on a µ which is still 5-15 times larger than the an-ticipated uncertainty from E989. Taken at face value, the A recent proposal for an independent determination of the HVPfrom the muon-electron elastic scattering has been put forwardin Ref. [4]. most recent lattice determination of the HVP [11] differsfrom the dispersive result by more than 3 standard devia-tions, and generates tensions with the global electroweakfits [12].All these facts call for an independent theoretically-sound lattice computation of the hadronic contributionto a µ at the per-mille level from first principles. Themain bottleneck toward this goal is the large statisticalerror in the Monte Carlo evaluation of the required corre-lation functions, see Ref. [3] and references therein. Theaim of this letter is to solve this problem by a novel com-putational paradigm based on multi-level Monte Carlointegration in the presence of fermions [13–15]. Withrespect to the standard approach, this strategy reducesthe variance exponentially with the temporal distance ofthe fields, thus opening the possibility of making negli-gible the contribution to the statistical error from long-distances. Here we focus on the HVP, but the strategyis general and can be applied to the HLbL, the isospin-breaking and electromagnetic contributions as well. II. THE SIGNAL-TO-NOISE PROBLEM
The HVP can be written as [16] a HVP µ = (cid:16) απ (cid:17) (cid:90) ∞ dx K ( x , m µ ) G ( x ) , (1)where α is the electromagnetic coupling constant, K ( x , m µ ) is a known analytic function which increasesquadratically at large x [17], m µ is the muon mass, and G ( x ) is the zero-momentum correlation function G ( x ) = (cid:90) d x (cid:104) J emk ( x ) J emk (0) (cid:105) (2) a r X i v : . [ h e p - l a t ] J u l of two electromagnetic currents J emk = i (cid:80) N f i =1 q i ¯ ψ i γ k ψ i .In this study we consider N f = 3, the three lighter quarksof QCD with the first two degenerate in mass, so that G ( x ) = G conn u,d ( x ) + G conn s ( x ) + G disc u,d,s ( x ) . (3)The light-connected Wick contraction G conn u,d ( x ) and thedisconnected one G disc u,d,s ( x ) are the most problematiccontributions with regard to the statistical error. In stan-dard Monte Carlo computations, the relative error of theformer at large time distances | x | goes as σ G connu , d ( x )[ G connu , d ( x )] ∝ n e M ρ − M π ) | x | , (4)where M ρ is the lightest asymptotic state in the iso-triplet vector channel, and n is the number of inde-pendent field configurations. For physical values of thequark masses, the difference ( M ρ − M π ) can be as largeas 3 . − . The computational effort, proportional to n , of reaching a given relative statistical error thus in-creases exponentially with the distance | x | . For the dis-connected contribution G disc u,d,s ( x ), the situation is evenworse since the variance is constant in time and thereforethe coefficient multiplying | x | is larger. At present thisexponential increase of the relative error is the barrierwhich prevents lattice theorists to reach a per-mille sta-tistical precision for the HVP. In order to mitigate thisproblem, in state-of-the-art calculations the contributionto the integral in Eq. (1) is often computed from MonteCarlo data only up to time-distances of 1 . − G ( x ), seeRef. [18] for an up-to-date review. III. MULTI-LEVEL MONTE CARLO
Thanks to the conceptual, algorithmic and technicalprogress over the last few years, it is now possible to carryout multi-level Monte Carlo simulations in the presenceof fermions [13, 14]. The first step in this approach is thedecomposition of the lattice in two overlapping domainsΩ and Ω , see e.g. Fig. 1, which share a common regionΛ . The latter is chosen so that the minimum distancebetween the points belonging to the inner domains Λ and Λ remains finite and positive in the continuum limit.Next step consists in rewriting the determinant of theHermitean massive Wilson-Dirac operator Q = γ D asdet Q = det (1 − w )det Q Λ det Q − det Q − , (5)where Q Λ , Q Ω , and Q Ω indicate the very same opera-tor restricted to the domains specified by the subscript.They are obtained from Q by imposing Dirichlet bound-ary conditions on the external boundaries of each do-main. The matrix w is built out of Q Ω , Q Ω and thehopping terms of the operator Q across the boundaries Λ Λ Λ Λ J emk J emk Ω Ω time s p a c e FIG. 1. Domain decomposition of the lattice adopted in thispaper. Periodic (anti-periodic) boundary conditions in thetime direction are enforced for gluons (fermions). in between the inner domains Λ and Λ and the commonregion Λ [14]. The denominator in Eq. (5) has alreadya factorized dependence on the gauge field since det Q Λ ,det Q − and det Q − depend only on the gauge field inΛ , Ω and Ω respectively.In the last step, the numerator in Eq. (5) is rewrittenas det (1 − w ) = det [1 − R N +1 (1 − w )] C (cid:81) N/ k =1 det (cid:2) ( u k − w ) † ( u k − w ) (cid:3) , (6)where u k and u ∗ k are the N roots of a polynomial approx-imant for (1 − w ) − , the numerator is the remainder, and C is an irrelevant constant. The denominator in Eq. (6)can be represented by an integral over a set of N/ and Λ [13–15] in-herited from w . When the polynomial approximation isproperly chosen, see below, the remainder in the numera-tor of Eq. (6) has mild fluctuations in the gauge field, andis included in the observable in the form of a reweightingfactor in order to obtain unbiased estimates.A simple implementation of these ideas is to dividethe lattice as shown in Fig. 1, where Λ and Λ have theshape of thick time-slices while Λ includes the remainingparts of the lattice. The short-distance suppression ofthe quark propagator implies that a thickness of 0 . is good enough,and is not expected to vary significantly with the quarkmass. This is the domain decomposition that we use forthe numerical computations presented in this letter.The Monte Carlo simulation is then performed usinga two-level scheme. We first generate n level-0 gaugefield configurations by updating the field over the entirelattice; then, starting from each level-0 configuration, wekeep fixed the gauge field in the overlapping region Λ ,and generate n level-1 configurations by updating thefield in Λ and in Λ independently thanks to the factor-ization of the action. The resulting gauge fields are thencombined to obtain effectively n · n configurations atthe cost of generating n · n gauge fields over the entirelattice. In particular, for each level-0 configuration, wecompute the statistical estimators by averaging the val-ues of the correlators over the n level-1 gauge fields. Pre-vious experience on two-level integration [13, 14, 21, 22] . . . . . . { σ / σ H M C } ( x , y ) x − y (fm) G conn u,d n = 1 n = 3 n = 3 ∗ n = 10 10 . . . a σ ( x , y ) x − y (fm) (cid:0) απ (cid:1) { KG conn u,d } × n = 1 n = 3 n = 3 ∗ n = 10 FIG. 2. Left: variance of the light-connected contraction as a function of the difference between the time-coordinates of thecurrents for n = 1 , ,
10 ( n = 3 ∗ is obtained by skipping 48 MDUs between two consecutive level-1 configurations). Dataare normalized to the analogous ones computed on CLS configurations generated by one-level HMC. Dashed lines representthe maximum reduction which can be obtained by two-level integration, namely 1 /n , in the absence of correlations betweenlevel-1 configurations. Grey bands indicate the thick time-slices where the gauge field is kept fixed during level-1 updates.Right: variance of the light-connected contribution to the integrand in Eq. (1). suggests that, with two independently updated regions,the variance decreases proportionally to 1 /n until thestandard deviation of the estimator is comparable withthe signal, i.e. until the level-1 integration has solvedthe signal-to-noise problem. From Eq. (4) we thus inferthat the variance reduction due to level-1 integration isexpected to grow exponentially with the time-distanceof the currents in Eq. (2). The overhead for simulat-ing the extra multi-boson fields increases the cost by anoverall constant factor which is quickly amortized by theimproved scaling. IV. LATTICE COMPUTATION
In order to assess the potential of two-level MonteCarlo integration, we simulate QCD with two dynamicalflavours supplemented by a valence strange quark. Glu-ons are discretized by the Wilson action while quarks bythe O ( a )-improved Wilson–Dirac operator, see Refs. [15,23] for unexplained definitions. Periodic and anti-periodic boundary conditions are imposed on the gluonand fermion fields in the time direction respectively, whileperiodic conditions are chosen for all fields in the spatialdirections. We simulate a lattice of size 96 × withan inverse bare coupling constant β = 6 /g = 5 .
3, corre-sponding to a spacing of a = 0 .
065 fm [23, 24]. The sizeof the lattice, rather large for a proof of concept, is cho-sen so to be able to accommodate a light pion and still bein the large volume regime, namely M π = 270 MeV and M π L ≥
4. The domains Λ and Λ are the union of 40consecutive time-slices, while each thick time-slice form-ing the overlapping region Λ is made of 8 time-slices cor-responding to a thickness of approximately 0 . N = 12.The very same action and set of auxiliary fields are usedeither at level-0 or level-1. The reweighting factor is es-timated stochastically with 2 random sources, which are enough for its contribution to the statistical error to benegligible. Further details on the algorithm and its im-plementation can be found in Ref. [15].We generate n = 25 level-0 configurations sepa-rated by 48 molecular dynamics units (MDU), so thatin practice they can be considered statistically uncorre-lated [23, 24]. For each level-0 background gauge field,we generate n = 10 configurations in Λ and Λ spacedby 16 MDUs. The connected contraction is calculatedby inverting the Wilson-Dirac operator on local sources,while the disconnected one is computed via split-evenrandom-noise estimators [25]. For each level-0 configu-ration, the statistical estimators are computed by aver-aging the correlators over the n combinations of level-1fields. The error analysis then proceeds as usual. Forthe sake of the presentation, we show results in physicalunits and properly renormalized: the central value of thelattice spacing is taken from Ref. [23], and the one of thevector-current renormalization constant from Ref. [26].We do not take into account their contributions to theerrors since on one side we are interested in investigatingthe statistical precision of the vector correlator computedvia two-level integration only, and on the other side thenumerical accuracy of those quantities can be improvedindependently.To single out the reduction of the variance due only totwo-level averaging, we carry out a dedicated calculationof correlation functions. We compute the light-connectedcontraction by averaging over 216 local sources put onthe time-slice ( y /a = 32) of Λ at a distance of 8 latticespacings from its right boundary and, as usual, by sum-ming over the sink space-position. This large number ofsources guarantees that the dependence of the varianceon the gauge field in the domains Λ and Λ is on equalfooting, since no further significant variance reductionis observed by increasing their number. We determinethe disconnected contraction by averaging each single-propagator trace over a large number of Gaussian randomsources, namely 768, so to have a negligible random-noisecontribution to the variance [25]. . . . x (fm) (cid:0) απ (cid:1) { KG } ( x ) × (fm − ) G conn u,d G conn s G disc u,d,s × − . . . x max0 (fm) a HVP µ ( x max0 ) × G conn u,d G conn s G disc u,d,s × FIG. 3. Left: best results for the contribution to the integrand in Eq. (1) from the light-connected (red squares), strange-connected (blue circles) and disconnected (green triangles) contractions as a function of the time coordinate. Right: bestresults for the contributions to a HVP µ from light-connected (red squares), strange-connected (blue circles), and disconnected(green triangles) contractions as a function of the upper extrema of integration x max0 . The variance of the light-connected contribution as afunction of the distance from the source is shown onthe left plot of Fig. 2. For better readability only thetime-slices belonging to Ω are shown, i.e. those rele-vant for studying the effect of two-level integration giventhe source position. Data are normalized to the varianceobtained with the same number of sources on CLS con-figurations which were generated with a conventionalone-level HMC [24, 27, 28]. The exponential reductionof the variance with the distance from the source is man-ifest in the data, with the maximum gain reached from2 . n = 10. The loss of about a factorbetween 2 and 3 with respect to the best possible scal-ing, namely n , either for n = 3 or 10 (dashed lines)is compatible with the presence of a residual correlationamong level-1 configurations. Indeed the variance reduc-tion for n = 3, obtained by skipping 48 MDUs betweenconsecutive level-1 configurations (labeled by n = 3 ∗ ), iscompatible with the n scaling at large distances withinerrors. In our particular setup, even for n = 10 the sta-tistical error at large distances scales de-facto with theinverse of the cost of the simulation rather than with itssquared root. This is easily seen by comparing the vari-ance reduction shown in the left plot of Fig. 2 with thecost of the simulation for n = 10. The latter is in fact 4times the one for n = 1 due to the different separationin MDU units between two consecutive configurations atlevel-0 and level-1.The power of the two-level integration can be betterappreciated from the right plot of Fig. 2, where we showthe variance of the light-connected contribution to theintegrand in Eq. (1) as a function of the time-distanceof the currents. The sharp rising of the variance com-puted by one-level Monte Carlo ( n = 1, red squares)is automatically flattened out by the two-level multi-boson domain-decomposed HMC ( n = 10, blue trian-gles) without the need for modeling the long-distancebehaviour of G conn u,d ( x ). https://wiki-zeuthen.desy.de/CLS/CLS. To further appreciate the effect of the two-level integra-tion, we compute the integral in Eq. (1) as a function ofthe upper extrema of integration x max0 which we allow tovary. For n = 1, the integral reads 446(26) and 424(38)for x max0 = 2 . . n = 10the analogous values are 467 . .
4) and 473 . . . . . n arerequired to render the variance approximately constant. V. RESULTS AND DISCUSSION
Our best result for the light-connected contributionto the integrand in Eq. (1) is shown on the left plotof Fig. 3 (red squares). It is obtained by a weightedaverage of the above discussed correlation functioncomputed on 32 point sources per time-slice on 7time-slices at y /a = { , , , , , , } and on216 sources at y /a = 32. We obtain a good statisticalsignal up to the maximum distance of 3 fm or so.The strange-connected contraction G conn s ( x ) is muchless noisy, and is determined by averaging on 16 pointsources at y /a = 32. Its value, shown on the left plot ofFig. 3 (blue circles), is at most one order of magnitudesmaller than the light contribution, and has a negligiblestatistical error with respect to the light one. Thebest result for the disconnected contribution has beencomputed as discussed in the previous section, and it isshown in the left plot of Fig. 3 as well (green triangles).It reaches a negative peak at about 1 . . a HVP µ · as a function of the upper extrema of inte-gration x max0 in Eq. (1). The light-connected part startsto flatten out at x max0 ∼ . x max0 = 3 . . . . x max0 = 3 . x max0 ∼ . − . x max0 = 3 . . x max0 = 3 . x max0 = 2 . a HVP µ = 522 . . · − .In this proof of concept study we have achieved a 1%statistical precision with just n · n = 250 configurationson a realistic lattice. This shows that for this light-quarkmass a per-mille statistical precision on a HVP µ is reach-able with multi-level integration by increasing n and n by a factor of about 4–6 and 2–4 respectively. Whenthe up and the down quarks becomes lighter, the gaindue to the multi-level integration is expected to increaseexponentially in the quark mass, hence improving evenmore dramatically the scaling of the simulation costwith respect to a standard one-level Monte Carlo. The change of computational paradigm presented here thusremoves the main barrier for making affordable, on thecomputer installations available today, the goal of aper-mille precision on a HVP µ .Here we focused on the main bottleneck in thecomputation of the HVP. It goes without saying thatthe very same variance-reduction pattern is expected towork out also for the calibration of the lattice spacing,the calculation of the electromagnetic corrections andthe HLbL.It is also interesting to notice that multi-level integra-tion can be well integrated with master-field simulationtechniques [29] if very large volumes turn out to benecessary to pin down finite-size effects at the per-millelevel. As a final remark, we stress that the very sameapproach is applicable to many other computationswhich suffer from signal-to-noise ratio problems, wherea similar breakthrough is expected [30]. ACKNOWLEDGMENTS
The generation of the configurations and the mea-surement of the correlators have been performedon the PC clusters Marconi at CINECA (CINECA-INFN, CINECA-Bicocca agreements, ISCRA B projectHP10BF2OQT) and at the Juelich Supercomputing Cen-tre, Germany (PRACE project n. 2019215140) whilethe R&D has been carried out on Wilson and Knuthat Milano-Bicocca. We thank these institutions andPRACE for the computer resources and the technicalsupport. We also acknowledge PRACE for awardingus access to MareNostrum at Barcelona SupercomputingCenter (BSC), Spain (n. 2018194651) where comparativeperformance tests of the code have been performed. Weacknowledge partial support by the INFN project “Highperformance data network”. [1] G. W. Bennett et al. (Muon g-2), Final Report of theMuon E821 Anomalous Magnetic Moment Measurementat BNL, Phys. Rev.
D73 , 072003 (2006), arXiv:hep-ex/0602035 [hep-ex].[2] J. Grange et al. (Muon g-2), Muon (g-2) Technical DesignReport, (2015), arXiv:1501.06858 [physics.ins-det].[3] T. Aoyama et al. , The anomalous magnetic mo-ment of the muon in the Standard Model, (2020),arXiv:2006.04822 [hep-ph].[4] G. Abbiendi et al. , Measuring the leading hadronic con-tribution to the muon g-2 via µe scattering, Eur. Phys.J. C , 139 (2017), arXiv:1609.08987 [hep-ex].[5] T. Blum, P. A. Boyle, V. G¨ulpers, T. Izubuchi, L. Jin,C. Jung, A. J¨uttner, C. Lehner, A. Portelli, and J. T.Tsang (RBC, UKQCD), Calculation of the hadronic vac-uum polarization contribution to the muon anomalousmagnetic moment, Phys. Rev. Lett. , 022003 (2018),arXiv:1801.07224 [hep-lat]. [6] D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo, andS. Simula, Electromagnetic and strong isospin-breakingcorrections to the muon g − , 114502 (2019), arXiv:1901.10462 [hep-lat].[7] C. Davies et al. (Fermilab Lattice, LATTICE-HPQCD,MILC), Hadronic-vacuum-polarization contribution tothe muon’s anomalous magnetic moment from four-flavor lattice QCD, Phys. Rev. D , 034512 (2020),arXiv:1902.04223 [hep-lat].[8] E. Shintani and Y. Kuramashi (PACS), Hadronic vac-uum polarization contribution to the muon g − latticeat the physical point, Phys. Rev. D100 , 034517 (2019),arXiv:1902.00885 [hep-lat].[9] A. G´erardin, M. C`e, G. von Hippel, B. H¨orz, H. B.Meyer, D. Mohler, K. Ottnad, J. Wilhelm, and H. Wit-tig, The leading hadronic contribution to ( g − µ fromlattice QCD with N f = 2 + 1 flavours of O( a ) im- proved Wilson quarks, Phys. Rev. D100 , 014510 (2019),arXiv:1904.03120 [hep-lat].[10] C. Aubin, T. Blum, C. Tu, M. Golterman, C. Jung, andS. Peris, Light quark vacuum polarization at the physicalpoint and contribution to the muon g −
2, Phys. Rev.
D101 , 014503 (2020), arXiv:1905.09307 [hep-lat].[11] S. Borsanyi et al. , Leading-order hadronic vacuum polar-ization contribution to the muon magnetic momentfromlattice QCD, (2020), arXiv:2002.12347 [hep-lat].[12] A. Crivellin, M. Hoferichter, C. A. Manzari, and M. Mon-tull, Hadronic vacuum polarization: ( g − µ versus globalelectroweak fits, (2020), arXiv:2003.04886 [hep-ph].[13] M. C`e, L. Giusti, and S. Schaefer, Domain decompo-sition, multi-level integration and exponential noise re-duction in lattice QCD, Phys. Rev. D93 , 094507 (2016),arXiv:1601.04587 [hep-lat].[14] M. C`e, L. Giusti, and S. Schaefer, A local factorizationof the fermion determinant in lattice QCD, Phys. Rev.
D95 , 034503 (2017), arXiv:1609.02419 [hep-lat].[15] M. Dalla Brida, L. Giusti, T. Harris, and M. Pepe, inpreparation, .[16] D. Bernecker and H. B. Meyer, Vector Correlators in Lat-tice QCD: Methods and applications, Eur. Phys. J.
A47 ,148 (2011), arXiv:1107.4388 [hep-lat].[17] M. Della Morte, A. Francis, V. G¨ulpers, G. Herdo´ıza,G. von Hippel, H. Horch, B. J¨ager, H. B. Meyer, A. Nyf-feler, and H. Wittig, The hadronic vacuum polarizationcontribution to the muon g − , 020, arXiv:1705.01775 [hep-lat].[18] V. G¨ulpers, Recent Developments of Muon g-2 from Lat-tice QCD, in (2019) arXiv:2001.11898 [hep-lat].[19] M. L¨uscher, A New approach to the problem of dynami-cal quarks in numerical simulations of lattice QCD, Nucl.Phys. B418 , 637 (1994), arXiv:hep-lat/9311007 [hep-lat].[20] A. Borici and P. de Forcrand, Systematic errors ofL¨uscher’s fermion method and its extensions, Nucl. Phys.
B454 , 645 (1995), arXiv:hep-lat/9505021 [hep-lat]. [21] M. L¨uscher and P. Weisz, Locality and exponential errorreduction in numerical lattice gauge theory, JHEP ,010, arXiv:hep-lat/0108014 [hep-lat].[22] M. Della Morte and L. Giusti, A novel approach for com-puting glueball masses and matrix elements in Yang-Millstheories on the lattice, JHEP , 056, arXiv:1012.2562[hep-lat].[23] G. P. Engel, L. Giusti, S. Lottini, and R. Sommer, Spec-tral density of the Dirac operator in two-flavor QCD,Phys. Rev. D91 , 054505 (2015), arXiv:1411.6386 [hep-lat].[24] P. Fritzsch, F. Knechtli, B. Leder, M. Marinkovic,S. Schaefer, R. Sommer, and F. Virotta, The strangequark mass and Lambda parameter of two flavor QCD,Nucl. Phys.
B865 , 397 (2012), arXiv:1205.5380 [hep-lat].[25] L. Giusti, T. Harris, A. Nada, and S. Schae-fer, Frequency-splitting estimators of single-propagatortraces, Eur. Phys. J.
C79 , 586 (2019), arXiv:1903.10447[hep-lat].[26] M. Dalla Brida, T. Korzec, S. Sint, and P. Vilaseca,High precision renormalization of the flavour non-singletNoether currents in lattice QCD with Wilson quarks,Eur. Phys. J.
C79 , 23 (2019), arXiv:1808.09236 [hep-lat].[27] L. Del Debbio, L. Giusti, M. L¨uscher, R. Petronzio, andN. Tantalo, QCD with light Wilson quarks on fine lattices(I): First experiences and physics results, JHEP , 056,arXiv:hep-lat/0610059 [hep-lat].[28] L. Del Debbio, L. Giusti, M. L¨uscher, R. Petronzio, andN. Tantalo, QCD with light Wilson quarks on fine lat-tices. II. DD-HMC simulations and data analysis, JHEP , 082, arXiv:hep-lat/0701009 [hep-lat].[29] A. Francis, P. Fritzsch, M. L¨uscher, and A. Rago, Master-field simulations of O(a)-improved lattice QCD: Algo-rithms, stability and exactness, Comput. Phys. Commun. , 107355 (2020), arXiv:1911.04533 [hep-lat].[30] L. Giusti, M. C`e, and S. Schaefer, Multi-boson block fac-torization of fermions, EPJ Web Conf.175