Pion Valence Quark Distribution from Matrix Element Calculated in Lattice QCD

TL;DR

This paper demonstrates a novel lattice QCD method to directly extract the pion valence quark distribution from coordinate space current-current correlations, offering a complementary approach to experimental global analyses.

Contribution

It introduces a new lattice QCD technique using antisymmetric vector and axial-vector currents to determine the pion valence quark distribution.

Findings

First exploratory lattice QCD calculation of pion valence distribution

Validation of the method with controllable systematics

Potential to complement experimental global analyses

Abstract

We present the first exploratory lattice QCD calculation of the pion valence quark distribution extracted from spatially separated current-current correlations in coordinate space. We show that an antisymmetric combination of vector and axial-vector currents provides direct information on the pion valence quark distribution. Using the collinear factorization approach, we calculate the perturbative tree-level kernel for this current combination and extract the pion valence distribution. The main goal of this article is to demonstrate the efficacy of this general lattice QCD approach in the reliable extraction of parton distributions. With controllable power corrections and a good understanding of the lattice systematics, this method has the potential to serve as a complementary to the many efforts to extract parton distributions in global analyses from experimentally measured cross sections. We perform our calculation on an ensemble of 2+1 flavor QCD using the isotropic-clover fermion action, with lattice dimensions $32^3\times 96$ at a lattice spacing \mbox{$a=0.127$ fm} and the quark mass equivalent to a pion mass $m_\pi \simeq 416$ MeV.

Tables

Table 1. Table 1: The lattice momenta p → = [ 0 , 0 , p z ] → 𝑝 0 0 subscript 𝑝 𝑧 \vec{p}=\left[0,0,p_{z}\right] of our interpolating operators and the momentum-smearing phases ζ 𝜁 \zeta applied for each lattice momenta in the direction d → = [ 0 , 0 , 1 ] → 𝑑 0 0 1 \vec{d}=\left[0,0,1\right] , as well as number of pion source points and source-sink separations. Quark sources comprising our interpolators were subsequently spatially smeared according to the Jacobi-smearing procedure with smearing parameters σ = 4.0 𝜎 4.0 \sigma=4.0 and n σ = 50 subscript 𝑛 𝜎 50 n_{\sigma}=50 .

$\overrightarrow{p}=[0,0,p_z]$	$\zeta$	No. of source points $(x_0,t)$	No. of source-sink    separations
$p=0.610\text{ GeV}$	$1.75$	$2$	$9$
$p=0.915\text{ GeV}$	$2.50$	$5$	$9$
$p=1.220\text{ GeV}$	$3.75$	$6$	$9$
$p=1.525\text{ GeV}$	$4.50$	$7$	$7$

Table 2. Table 2: The fit parameters of the V-A matrix elements with ξ → = [ 0 , 0 , a ] → 𝜉 0 0 𝑎 \vec{\xi}=\left[0,0,a\right] for momenta along the z 𝑧 z -direction p = 0.610 𝑝 0.610 p=0.610 GeV and p = 1.525 𝑝 1.525 p=1.525 GeV, respectively.

$p$ [GeV]	$\xi$	$A$	$B$	${\Delta }_{\mathrm{eff}}$	${\chi }^2$/d.o.f.
0.610 GeV	$1a$	$0.102\left(5\right)$	-0.028(11)	0.054(20)	1.21
0.610 GeV	$4a$	$0.083\left(4\right)$	-0.026(13)	0.062(42)	0.89
1.525 GeV	$1a$	$0.097\left(8\right)$	-0.267(513)	0.809(503)	0.15
1.525 GeV	$4a$	$0.049\left(9\right)$	-0.107(148)	0.450(420)	0.19

Equations

$${\sigma }_n(\xi ,p)=⟨p|T\{{\mathcal{O}}_n(\xi )\}|p⟩,$$

$${\mathcal{O}}_{j_1{j}_2}(\xi )\equiv$$

$${\sigma }_{j_1{j}_2}(\omega ,{\xi }^2)=⟨p|T\{{\mathcal{O}}_{j_1{j}_2}(\xi )\}|p⟩,$$

$${\mathcal{O}}_{VV}^{\mu \nu }(\xi )\equiv$$

$${\sigma }_{VV}^{\mu \nu }(\xi ,p)=p^{\mu }{p}^{\nu }{T}_1(\omega ,{\xi }^2)+\frac{1}{2}(p^{\mu }{\xi }^{\nu }+{\xi }^{\mu }{p}^{\nu })T_2(\omega ,{\xi }^2)$$

$$+g^{\mu \nu }{T}_3(\omega ,{\xi }^2)+{\xi }^{\mu }{\xi }^{\nu }{T}_4(\omega ,{\xi }^2)$$

$${\tilde{\sigma }}_n(\tilde{\omega },q^2)\equiv \int \frac{d^4\xi }{{\xi }^4}{e}^{iq\cdot \xi }{\sigma }_n(\omega ,{\xi }^2),$$

$${\sigma }_n^h(\omega ,{\xi }^2,p^2)$$

$${\sigma }_n^{q\left(0\right)}\left(\omega ,{\xi }^2\right)$$

$${\sigma }_n^{q\left(0\right)}\left(\omega ,{\xi }^2\right)=K_n^{q\left(0\right)}\left(\omega ,{\xi }^2\right).$$

$${\mathcal{O}}_{ij}^{\mu \nu }\left(\xi \right)={\xi }^4{\mathcal{J}}_i^{\mu }\left(\xi /2\right){\mathcal{J}}_j^{\nu }\left(-\xi /2\right),$$

$${\sigma }_{ij}^{\mu \nu }\left(\xi ,p\right)$$

$${\mathcal{M}}_{ij}^{(a)}$$

$${\mathcal{M}}_{ji}^{(b)}=\frac{{\xi }^4}{2}{e}^{-ik\cdot \xi }\mathrm{Tr}\left[\left(\gamma \cdot k\right){\Gamma }_j^{\nu }\int \frac{d^{4}l}{{\left(2\pi \right)}^4}\frac{-i\gamma \cdot l}{l^2+iϵ}{e}^{-il\cdot \xi }{\Gamma }_i^{\mu }\right]$$

$${\sigma }_{ij}^{\mu \nu \left(0\right)}\left(p\cdot \xi ,p;x,\xi \right)$$

$${\sigma }_{ij}^{\mu \nu }\left(\xi ,p\right)=⟨\pi \left(p\right)|\left(\mathcal{P}\mathcal{T}\right){\left({\mathcal{O}}_{ij}^{\mu \nu }\left(\xi \right)\right)}^{†}{\left(\mathcal{P}\mathcal{T}\right)}^{-1}|\pi \left(p\right)⟩.$$

$$\left(\mathcal{P}\mathcal{T}\right){\mathcal{J}}_A^{\mu }\left(\xi \right){\left(\mathcal{P}\mathcal{T}\right)}^{-1}=-{\mathcal{J}}_A^{\mu }\left(-\xi \right)$$

$$\left(\mathcal{P}\mathcal{T}\right){\mathcal{J}}_V^{\mu }\left(\xi \right){\left(\mathcal{P}\mathcal{T}\right)}^{-1}={\mathcal{J}}_V^{\mu }\left(-\xi \right)$$

$$\frac{1}{2}\left[{\sigma }_{VA}^{\mu \nu }\left(\xi ,p\right)+{\sigma }_{AV}^{\mu \nu }\left(\xi ,p\right)\right]$$

$$\equiv {ϵ}^{\mu \nu \alpha \beta }{\xi }_{\alpha }{p}_{\beta }{T}_1\left(\omega ,{\xi }^2\right)+\left(p^{\mu }{\xi }^{\nu }-{\xi }^{\mu }{p}^{\nu }\right)T_2\left(\omega ,{\xi }^2\right)$$

$$T_1\left(\omega ,{\xi }^2\right)$$

$$T_2\left(\omega ,{\xi }^2\right)$$

$$T_1\left(\omega ,{\xi }^2\right)$$

$$T_2\left(\omega ,{\xi }^2\right)$$

$$T_i\left(\omega ,{\xi }^2\right)$$

$${\sigma }_{VA}^{\mu \nu \left(0\right)}$$

$$T_1^{q(0)}(x\omega ,{\xi }^2)$$

$$T_2^{q(0)}(x\omega ,{\xi }^2)$$

$$C_1^{q(0)}(x\omega ,{\xi }^2)$$

$$C_2^{q(0)}(x\omega ,{\xi }^2)$$

$${\tilde{T}}_1$$

$$C_{2pt}(p,T)=⟨{\Pi }_p(T){\overline{\Pi }}_p(0)⟩,$$

$$C_{2pt}(p,T)=\sum _n\frac{|Z_n{|}^2}{2E_n(p)}{e}^{-E_n(p)T},$$

Full text

††preprint: JLAB-THY-19-2847

I Introduction

In the hard scattering processes involving hadrons, such as in the deep inelastic scattering (DIS) of leptons on hadrons, the experimentally measured cross sections are a combination of short- and long-distance physics. The inclusive DIS cross section can be factorized into a short-distance partonic hard part which is calculable order by order in perturbation theory and a long-distance hadronic part which can be represented by universal and nonperturbative distribution functions, called the parton distribution functions (PDFs), plus corrections suppressed by inverse power of large momentum transfer of the scattering. It is the QCD factorization theorem Collins:1989gx which enables us to connect the dynamics of quarks and gluons to the physically measured hard scattering cross sections of identified hadrons. The collinear (CO) divergences of the partonic scattering are absorbed into the nonperturbative PDFs, leaving an infrared-safe and perturbatively computable hard contribution. According to Feynman’s parton model Feynman , the unpolarized PDFs give the probability to find partons [i.e. quark ($q$), antiquark ($\bar{q}$), gluon ($g$)] in a hadron as a function of the fraction $x$ of the hadron’s longitudinal momentum carried by the parton, probed at a factorization scale $\mu$. For example, if a parton of type $i$ carries a fraction $x$ of hadron’s momentum, then the probability to find the parton is given by $f_{i}(x,\mu^{2})dx$. An accurate and precise knowledge of parton distribution functions is required for the cross section predictions of both Standard Model and Beyond Standard Model processes at existing and future particle colliders, such as the LHC and Electron-Ion Collider (EIC). The precision of numerous experimentally measured observables, such as the $W$-boson mass, weak-mixing angle and Higgs cross section, is driven by a detailed knowledge and precision of PDFs.

Since a precise knowledge of PDFs is required for the analysis and interpretation of scattering experiments, as discussed above, considerable effort has been made to determine PDFs and their uncertainties by global fitting collaborations such as MMHT Harland-Lang:2014zoa , CT Dulat:2015mca , NNPDF Ball:2017nwa , HERAPDF Alekhin:2017kpj , and JAM Ethier:2017zbq . PDFs, the catalyst of many observables in hadronic scattering, are becoming better determined as experimental data sets increase and the global analysis community implements more sophisticated schemes to quantify systematic uncertainties.

The valence quark distribution of the pion is of particular theoretical interest, as the pion is the lightest QCD bound state and the Goldstone mode associated with dynamical chiral symmetry breaking. The pion PDF has been measured through pionic Drell-Yan experiments at CERN Badier:1983mj ; Betev:1985pf and Fermilab Conway:1989fs . Several analyses in Refs. Owens:1984zj ; Aurenche:1989sx ; Sutton:1991ay ; Gluck:1991ey ; Wijesooriya:2005ir ; Aicher:2010cb ; Barry:2018ort of these experimental data have been performed to determine the pion valence distribution. Among these analyses, it has been emphasized in Ref. Aicher:2010cb that the next-to-leading-logarithmic threshold resummation effects in the calculation of the Drell-Yan cross section are important and give a softer valence distribution which falls of as $(1-x)^{2}$ near $x\to 1$, consistent with the prediction based on the framework of perturbative QCD in Refs. Farrar:1979aw ; Berger:1979du ; Brodsky:1994kg . There are also different model predictions for the large-$x$ behavior of the pion valence distribution, some of which predict a harder fall-off as $(1-x)$ Shigetani:1993dx ; Davidson:1994uv ; Melnitchouk:2002gh ; deTeramond:2018ecg or $(1-x)^{2}$, such as in Dyson-Schwinger type models Hecht:2000xa ; Chen:2016sno . Therefore, an ab initio knowledge of the correct large-$x$ behavior of the pion valence PDF can serve as a discriminator of different model calculations. In this paper, we will present a calculation of the pion valence PDF using “lattice cross sections” proposed in Refs. Ma:2014jla ; Ma:2017pxb . In this approach, one factorizes a hadronic matrix element, such as a two current correlator, into the PDFs and short-distance matching coefficients, from which PDFs could be extracted from lattice calculated hadronic matrix element with the perturbatively calculated matching coefficients.

Lattice QCD is a Monte Carlo method for numerically evaluating QCD in a finite, discretized Euclidean spacetime. To date, Lattice QCD has emerged as the most rigorous and systematic tool for studying QCD nonperturbatively. As introduced by Feynman Feynman:1969ej , PDFs are defined through light cone matrix elements of certain bilocal operators. These light cone matrix elements cannot be directly calculated on the Euclidean lattice because the light cone collapses to a point in Euclidean spacetime. Recently several methods have been introduced to go beyond the calculations of the first few moments of PDFs on the lattice, such as the path-integral formulation of the deep-inelastic scattering hadronic tensor Liu:1993cv ; Liu:1999ak , the inversion method Horsley:2012pz , quasi-PDFs Ji:2013dva , and pseudo-PDFs Radyushkin:2017cyf to obtain the $x$-dependent hadron structure functions. A coordinate-space method for the calculation of light-cone distribution amplitudes has also been employed Braun:2007wv . Significant achievements in the lattice QCD implementations of these approaches have been made in recent years Bali:2018spj ; Chambers:2017dov ; Alexandrou:2018pbm ; Lin:2018qky ; Orginos:2017kos . References to many other lattice QCD calculations, the current status and challenges for a meaningful comparison of these lattice calculations with the global fits of PDFs can be found in Refs. Lin:2017snn ; Monahan:2018euv ; Cichy:2018mum . The readers are also referred to several lattice QCD calculations of lower moments of pion PDF in Refs. Best:1997qp ; Detmold:2003tm ; Guagnelli:2004ga ; Capitani:2005jp ; Bali:2013gya ; Abdel-Rehim:2015owa ; Oehm:2018jvm . Recently, the relation of moments of Ioffe time parton distribution functions to matrix elements of nonlocal operators computed in lattice QCD has been studied in Karpie:2018zaz .

The remainder of this article is organized as follows.

In Secs. II and III we discuss what are the good lattice cross sections, and the essence of calculation of these lattice cross sections in coordinate space. In Sec. IV, we present the derivation of the tree-level perturbative kernel for an antisymmetric vector and axial-vector current combination, and show how one can factorize the associated hadronic matrix element to extract pion valence distribution in a lattice QCD calculation. We present the numerical methods and results in Secs. V and VI. We compare the pion valence distribution extracted in this calculation with other calculations and different fits of the experimental data. Finally, we summarize our results and outline the future directions of this method to obtain pion valence quark distribution with controlled systematics.

II “Good” Lattice Cross Sections

In similarity with the extraction of PDFs in a global fit through the factorization of different experimental cross sections, a method was proposed Ma:2014jla to extract PDFs from lattice QCD calculations of hadronic matrix elements, called lattice cross sections (LCSs), which in the framework of QCD can be factorized into a perturbatively calculable hard part and the nonperturbative PDFs with a small and controllable power correction. Analogous to the cross sections measured in an experiment, these hadronic matrix elements computed on the lattice are constructed to be time independent, defined by equal-time operators and have a well-defined continuum limit; hence the name lattice “cross sections.” It has been shown in Ref. Briceno:2017cpo that as long as such operators have no temporal extent, a matrix element calculated in Euclidean space will equal its counterpart in Minkowski space.

To be more specific, such hadronic matrix elements are good lattice cross sections if they have the following properties:

• •

are a Lorentz covariant single-hadron matrix element computable on a Euclidean lattice,

• •

have a well-defined continuum limit when the lattice spacing $a\to 0$,

• •

are factorizable to PDFs convoluted with infrared (IR)-safe hard coefficients, plus a small and controllable power correction.

The single-hadron matrix elements of renormalized nonlocal operators ${\cal O}_{n}({\xi})$ can be written as, suppressing renormalization scale dependence,

[TABLE]

where the subscript $n$ is a label for different operators, $T$ stands for time-ordering, $p$ the hadron momentum, and $\xi$ ($\xi^{2}\neq 0$ ) is the largest separation of all fields in the operator ${\cal O}_{n}$. One choice which allows for factorization is a pair of parton field operators linked by a Wilson line operator and has been used for the calculation of pseudo-PDFs Radyushkin:2017cyf and quasi-PDFs Ji:2013dva . A broader class of operators with factorizable matrix elements is pairs of spacelike separated currents.

The operator can be chosen to be a Lorentz scalar, such as

[TABLE]

where $j_{1}$ and $j_{2}$ are currents with no Lorentz indices such as $\bar{\psi}\psi$ or $\bar{\psi}\not{\xi}\psi$, $d_{j}$ and $Z_{j}$ are the dimension and renormalization constant of the current $j$, respectively, and the overall dimensional factor is introduced so that the matrix elements in Eq. (1) are dimensionless with our normalization, $\langle p|p^{\prime}\rangle=(2E_{p})(2\pi)^{3}\delta^{3}(p-p^{\prime})$. In this case the LCS can be written in terms of only Lorentz invariants

[TABLE]

where the Lorentz scalar $\omega\equiv p\cdot\xi$ is the Ioffe time Ioffe:1969kf , and $p^{2}$ dependence is suppressed. This type of LCS can be directly factorized into the PDF when the separation $|\xi|$ is small Ma:2014jla .

Other choices of operators can have a more complicated Lorentz structure, such as the vector-vector matrix element of the operator

[TABLE]

where $j_{V}^{\mu}$ is the vector current which requires no renormalization constants. This type of LCS will need to be decomposed into the Lorentz structures allowed by symmetry. The functions of Lorentz invariants which accompany these Lorentz structures are the objects which will be factorized into PDFs. For the case of the operator ${\cal O}^{\mu\nu}_{VV}$, its matrix element of an unpolarized hadron state is symmetric in $\{\mu,\nu\}$ and can be decomposed as,

[TABLE]

where $T_{i}$ are Lorentz invariant functions. In the following section, we discuss the vector and axial-vector current combination which will be used to extract the pion valence quark distribution.

It is worth noting that the LCSs have the following analogs to hard scattering experiments:

• •

the label “$n$” in Eq. (1) is related to the dynamical features of LCSs and mimics different processes in experiments.

• •

$p$ and $\xi$ are analogous to observed scales defining the collision kinematics; $p$ relating to the collision energy $\sqrt{S}$ and $\xi^{2}$ relating to the hard probe $\frac{1}{Q^{2}}$.

III The Essence of Calculation in Coordinate Space

It is crucial to mention that a large $p$ alone does not guarantee the applicability of the operator product expansion of the matrix element and contributions from large $\xi$ can invalidate the perturbative factorization, whether for the case of quark-antiquark fields linked by a Wilson line, or for the case of spatially-separated currents. The validity of perturbative factorization requires that the separation’s scale, $\xi^{2}$, be much smaller than the inverse square of typical hadronic scale, $\Lambda_{\rm QCD}$, namely $\xi^{2}\Lambda^{2}_{\rm QCD}\ll 1$.

We now show why the coordinate space approach provides distinct theoretical advantages. One can write the Fourier transform of $\sigma_{n}(\omega,\xi^{2})$:

[TABLE]

where corresponding ${\cal O}_{n}$ can be any operator defined in Eq. (1), and $\widetilde{\omega}\equiv\frac{2p\cdot q}{-q^{2}}=\frac{1}{x_{\rm B}}$ with $x_{\rm B}$ the Bjorken variable for the lepton-hadron DIS. However, though $q$ is related to $\xi$ through the Fourier transform above, it is not a one-to-one relation. $\tilde{\sigma}$ with small or large values of $q$ will receive contribution from $\sigma$ with all values of $\xi$, small and large. In particular, it can involve contributions from values of $\xi\gg\frac{1}{\Lambda_{\rm QCD}}$, thereby violating factorization. This is why, in the LCSs approach, $\xi$ is a very well-defined quantity, analogous to a hard probe of hadron structure in a DIS experiment. It has also been demonstrated in Ref. Ma:2017pxb that, the non-analytic cut of $\widetilde{\sigma}_{n}$ comes from the integration region of large $\xi$. That is, even if we demand $|q^{2}|\gg\Lambda_{\rm QCD}^{2}$, $\widetilde{\sigma}_{n}$ in momentum space can always receive contribution from large $\xi$ region so long as $\widetilde{\omega}^{2}>1$. On the other hand, in coordinate space, if we fix $\xi$ to be short-distance, we do not have contribution from the large $\xi$ region and thus ${\sigma}_{n}$ has a good analytic behavior.

IV Factorization

The lattice calculable hadronic matrix elements of Eq. (1) are shown in Ref. Ma:2017pxb to be factorizable into PDFs with perturbatively calculable coefficients by applying the operator product expansion (OPE) to the nonlocal operator $\mathcal{O}_{n}\left(\xi\right)$, with small but nonvanishing $\xi^{2}$

[TABLE]

where $\sigma_{n}^{h}$ is $\mathcal{O}_{n}\left(\xi\right)$ measured in a hadron $h$, $K_{n}^{a}$ are the parton flavor $a\in\{q,\overline{q},g\}$ contributions to the perturbative hard coefficients with corresponding PDF $f^{h}_{a}$, and factorization scale $\mu^{2}$. The short-distance coefficient functions $K_{n}^{a}\left(\omega,\xi^{2};p^{2},\mu^{2}\right)$ are determined by applying the factorized formula in Eq. (LABEL:eq:factor) to an asymptotic parton of momentum $p$ with $p^{2}=0$ and flavor $a=q,\overline{q}$, or $g$, and expanding each side as a power series in the strong coupling constant $\alpha_{s}$. Although the perturbative coefficient functions are process-dependent, they apply equally to different external hadron states. At leading-order $\mathcal{O}\left(\alpha_{s}\right)$, the matching coefficient only receives quark contributions and the factorized relation in Eq. (LABEL:eq:factor) becomes

[TABLE]

where $p^{2}=0$ is suppressed and $f_{a}^{q\left(0\right)}\left(x,\mu^{2}\right)=\delta\left(1-x\right)\delta^{qa}$ is the quark distribution of an asymptotic quark at zeroth order in $\alpha_{s}$ and does not have the factorization scale $\mu^{2}$-dependence. Upon substitution of $f_{a}^{q\left(0\right)}$ into Eq. (7), it can be shown

[TABLE]

where the active quark momentum $p^{2}=0$ is suppressed. Thus determination of the leading-order coefficient functions $K_{n}^{q\left(0\right)}$ follows directly from the expression of the matrix elements $\sigma_{n}^{q\left(0\right)}\left(\omega,\xi^{2}\right)$ in the coordinate space.

Specializing to the case of the pion, consider a generic tensor operator

[TABLE]

and have it evaluated in the pion state $\ket{\pi\left(p\right)}$, where $\mathcal{J}_{k}$ is a local quark bilinear and $\xi^{4}$ is included to maintain an overall dimensionless matrix element. By examining the path-integral definition of an arbitrary operator defined at a single Euclidean time removes complications in analytically continuing our results back to Minkowski space. In this case, the general time-ordering of $\mathcal{O}_{ij}^{\mu\nu}\left(\xi\right)$ is instead expressed as a sum of diagrams with momenta flowing in/out of the fixed current locations. We define the matrix element of $\mathcal{O}_{ij}^{\mu\nu}\left(\xi\right)$ in the pion as

[TABLE]

Projecting onto an asymptotic quark state, we are left with two distinct diagrams at leading order (LO):

Depending on the current-current combinations considered, the resulting Lorentz decomposition of $\sigma^{\mu\nu}_{ij}\left(\xi,p\right)$ will introduce numerous scalar form factors consistent with parity and time-reversal invariance. It is these form factors that will provide information on a wide array of distribution functions, when factorized according to Eq. (LABEL:eq:factor). A general expression from which the LO perturbative kernels can be obtained follows from application of perturbative formulae to the diagrams above. Averaging over quark spin, the ordering depicted in Fig. (1) yields

[TABLE]

where an inverse Fourier transform has been used to express the quark propagator from $-\xi/2\rightarrow\xi/2$ in coordinate space. The second ordering, shown in Fig. (1), similarly yields

[TABLE]

Combining Eqs. (11) and (12) and writing the quark momentum as $k_{\mu}=xp_{\mu}$, we obtain a general relation in the LO denoted by the superscript $(0)$ as

[TABLE]

from which the kernels $K_{n}^{q\left(0\right)}\left(\omega,\xi^{2};x\right)$ with $\omega=p\cdot\xi$ can be isolated for currents $\{i,j\}$.

Given invariance of the strong interaction under parity ($\mathcal{P}$) and time-reversal ($\mathcal{T}$) transformations, the pion matrix element $\sigma_{ij}^{\mu\nu}\left(\xi,p\right)$ has the following property,

[TABLE]

In this work, we consider the case of a vector $\mathcal{J}_{V}^{\mu}=\overline{\psi}\gamma^{\mu}\psi$ and axial-vector $\mathcal{J}_{A}^{\nu}=\overline{\psi}\gamma^{\nu}\gamma^{5}\psi$ current combination, whose transformation properties are

[TABLE]

With these transformation properties, we find that the following combination of these two currents, $\sigma_{VA}^{\mu\nu}\left(\xi,p\right)+\sigma^{\mu\nu}_{AV}\left(\xi,p\right)\equiv\bra{\pi\left(p\right)}\left[\mathcal{O}^{\mu\nu}_{VA}\left(\xi\right)+\mathcal{O}_{AV}^{\mu\nu}\left(\xi\right)\right]\ket{\pi\left(p\right)}$, is antisymmetric in Lorentz indices, $\{\mu,\nu\}$, and can be expressed in terms of two dimensionless scalar form factors as

[TABLE]

where $T_{i}\left(\omega,\xi^{2}\right)$ are the dimensionless functions of the Lorentz invariants $\{\omega,\xi^{2}\}$.

The dimensionless functions are isolated by taking appropriate tensor contractions of the antisymmetric matrix element in Eq. (15),

[TABLE]

A judicious choice of $\xi$, $p$, and Lorentz indices $\{\mu,\nu\}$, exposes the structure functions $T_{1}\left(\omega,\xi^{2}\right)$ and $T_{2}\left(\omega,\xi^{2}\right)$ without recourse to a full tensor contraction as in Eqs. (IV) and (IV). To isolate the structure functions we stipulate $p=(p^{0},0,0,p^{3})$ and $\xi=(0,0,0,\xi^{3})$. $T_{1}$ is then isolated by choosing $\mu=1$ and $\nu=2$:

[TABLE]

While $T_{2}$ is isolated by choosing $\mu=0$ and $\nu=3$:

[TABLE]

From Eq. (LABEL:eq:factor), the $T_{i}\left(\omega,\xi^{2}\right)$ can be thus factorized as

[TABLE]

with the perturbatively calculable matching coefficients $C^{a}_{i}\left(x\omega,\xi^{2},\mu^{2}\right)$ for parton flavor $a$.

Similar to the derivation of Eq. (13), the LO contribution to the antisymmetric pion matrix element in Eq. (15) is given by

[TABLE]

with {\rm Tr}{\big{(}\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}\gamma^{5}\big{)}}=-4i\epsilon^{\mu\nu\rho\sigma} dictated by the convention $\epsilon^{0123}=1$. Substituting Eq. (21) into Eqs. (IV) and (IV), we obtain the two scalar form factors of a quark state of momentum $k=xp$ at the LO, respectively,

[TABLE]

With $f_{a}^{q\left(0\right)}\left(x,\mu^{2}\right)=\delta\left(1-x\right)\delta^{qa}$, we obtain the lowest order perturbative coefficients in Eq. (20) as,

[TABLE]

respectively.

We have the LO momentum-space scalar form factors by performing a Fourier transformation in $\omega$,

[TABLE]

where $q(-\tilde{x})=-\overline{q}(\tilde{x})$ is used, $q_{\rm v}\!\left(\tilde{x}\right)\equiv[q\left(\tilde{x}\right)-\overline{q}\left(\tilde{x}\right)]$ is the valence quark distribution, and the $\xi^{2}$ or the factorization scale dependence is suppressed since we are working at the LO approximation. Equation (26) implies that $\tilde{T}_{1}(\tilde{x},\xi^{2})$ is proportional to the valence quark PDF with momentum fraction $\tilde{x}$, which is actually true to all orders due to the symmetry of the coefficient function $C_{1}^{q}(x\omega,\xi^{2},\mu^{2})=-C_{1}^{q}(-x\omega,\xi^{2},\mu^{2})$.

Therefore direct information on the pion’s valence quark distribution $q_{\rm v}\left(\tilde{x},\mu^{2}\right)$ is accessible by evaluating the antisymmetric combination of vector and axial-vector (V-A) current-current correlators, up to an overall factor of $1/\pi^{2}$ and corrections in powers of $\alpha_{s}$ and/or $\xi^{2}\Lambda_{\rm QCD}^{2}$.

It has been shown in Ref. Ma:2017pxb that the validity of operator product expansion (OPE) guarantees that $T_{1}$ is an analytic function of $\omega$, as is its Taylor series around $\omega=0$. By keeping $\xi$ to be short distance and increasing $\omega$ by increasing $p$, there exists no way for new divergences to appear in $T_{1}$. Therefore, $T_{1}$ remains an analytic function of $\omega$ unless $\omega=\infty$ and the factorization holds for any values of $\omega$ and $\xi^{2}$ as long as $\xi$ is short distance, similar to the scenario of the factorization of experimental cross sections.

V Numerical Methods

This calculation is performed on a lattice gauge ensemble of 490 configurations generated by the JLab/W&M Collaboration lattices . This ensemble employs 2+1 flavors of clover Wilson fermions and a tree-level tadpole improved Symanzik gauge action. The strange quark mass was set by requiring the ratio $\left(2M^{2}_{K^{+}}-M^{2}_{\pi^{+}}\right)/M_{\Omega^{-}}$ to assume its physical value. The configurations were generated using a rational Hybrid Monte Carlo update algorithm rhmc . The fermion action includes a single iteration of stout smearing with weight $\rho=0.125$. This smearing makes the employed tadpole corrected tree-level clover coefficient, $c_{sw}$, very close to the nonperturbative value determined, a posteriori, by the Schrödinger functional method.

The extraction of hadron-to-hadron matrix elements in lattice QCD requires the calculation of correlation functions. The 2-point function is a vacuum expectation value of two interpolating fields separated in Euclidean time $T$:

[TABLE]

where the interpolating field $\Pi_{p}$ is an operator with quantum numbers of a pion with momentum $p$. A spectral decomposition of the 2-point function is given by the following tower of exponentials

[TABLE]

where the sum is over all energy eigenstates $n$ with quantum numbers of the pion, $Z_{n}=\langle 0|\Pi_{p}|n\rangle$ is the overlap factor between the operator and the $n{\rm th}$ excited state and $E_{n}(p)$ is the energy of that state with momentum $p$. In the large Euclidean time limit, this correlation function will be dominated by the ground state.

A good choice of interpolating field will have a large overlap factor with the ground state while simultaneously having poor overlap with excited states. For low-momenta or states at rest, spatial smearing is a well-established method to reduce the overlap of pointlike interpolators onto high energy eigenstates. We employ in this work the Jacobi-smearing procedure Allton:1993wc , in which pointlike quark fields are smeared according to

[TABLE]

where $\nabla^{2}$ is the three dimensional gauge-covariant discretization of the Laplacian, $\sigma$ the smearing “width” and $n_{\sigma}$ the number of applications of the smearing kernel onto the pointlike quark fields. For highly-boosted states, however, the overlap of even spatially-smeared interpolators can become suboptimal. To ameliorate the effects of excited-states and improve the overlap of our interpolators onto boosted pions, we implement a combination of the Jacobi and momentum-smearing Bali:2016lva techniques. In practice we apply appropriately constructed phases to the underlying gauge fields prior to source creation according to

[TABLE]

where $\vec{d}$ is the direction in which phases are applied, with magnitude $\zeta$ tuned for each desired momenta. The final interpolating fields are given by

[TABLE]

where $\tilde{q}$ is a light quark field constructed with the combined application of momentum smearing and Jacobi smearing. The smearing parameters used were varied for each momentum and are shown in Table I. Due to the decreasing signal to noise ratio, the higher momentum states required more source points and shorter time separation between the pion operators. These values are also shown in Table I.

For the calculation of any good LCS, the composite operators used have finite spatial extent $\xi$. Introduction of a heavy auxiliary quark field $Q$ ($m_{Q}>m_{l}$), such that our operators are of the form $\mathcal{O}(t)=\mathcal{J}_{\Gamma}^{\dagger}(\xi,t)\mathcal{J}_{\Gamma^{\prime}}(0,t)$ with $\mathcal{J}_{\Gamma}=\bar{q}\Gamma Q$, limits the available phase space between the two currents thereby reducing the statistical noise. An auxiliary heavy quark has also been used in Ref. Detmold:2005gg to remove higher twist contamination in the calculation of moments of the PDF and the distribution amplitude (DA). For our calculation of the pion valence distribution, multiple auxiliary quark masses between the light and strange quark mass were tested. A slight improvement in the signal-to-noise ratio from the heavier masses was observed for the larger momenta. We set the auxiliary quark propagator to the strange quark mass for the remainder of this calculation. In addition, to minimize excited state contamination, the operator insertion time ($t$) will be fixed to be midway between the source and sink interpolators (i.e. $t=\frac{T}{2}$).

The 4-point correlation function is constructed using a modified sequential source technique. Because we are not performing a time slice momentum projection at the operator, the standard sequential source method using the operator as sequential source does not work here. However, for the case of meson there is a straightforward implementation where momentum projections are performed at source and sink meson operators and the corresponding correlation functions are computed as chain of sequential sources as described below. The correlation function is expressed as follows

[TABLE]

where $(x_{0},t)$ is a randomly determined source point, $G_{Q}(y^{\prime};x^{\prime})$ is the flavor $Q$ auxiliary quark propagator from $x^{\prime}$ to $y^{\prime}$, and $I^{p}_{q}(y^{\prime};x^{\prime})$ is the modified sequential source with flavor $q$-quarks and pions at momentum $p$. The modified sequential source is constructed through sequential inversions of the light quark Dirac operator, reusing already calculated propagators. Heuristically, the modified sequential source is constructed by calculating the light quark propagator from a point-source located at one of the currents to the source interpolator, using this object as a source for a subsequent propagator to the sink interpolator, and lastly using this larger object as a source for propagation from the sink to the second current. This construction is done by solving the following sequence of systems of equations for $G_{q}$, $H^{p}_{q}$, and $I^{p}_{q}$.

[TABLE]

where $D_{q}$ is the Dirac matrix for the light quarks and $S(x;x^{\prime})$ represents the smearing procedure. The phase projects the interpolating fields onto definite momentum, while the signs of the momentum-smearing phases applied to the quark fields must be treated carefully to ensure $\Pi$ and $\overline{\Pi}$ correctly project onto states with a given momentum. Note that the smearing procedures are applied once for each of the quark fields in the interpolating fields, and no smearing is applied at the current insertions. This procedure is shown diagrammatically in Fig 2. An advantage of this procedure is that the correlation function for any current pair with any separation can be calculated without additional Dirac matrix inversions. On the other hand, for each choice of momentum and source-sink separations, two additional light quark Dirac matrix inversions are required per gauge configuration, denoted by
1 and
2 in Eq. (V).

VI Numerical Results and Extraction of Pion Valence Distribution

Reliable extraction of hadronic matrix elements, in part, hinges on how well a lattice calculation can systematically quantify and reduce excited-state contamination. In this section, we present the numerical results of our calculation of matrix elements of spatially-separated antisymmetric V-A currents. As discussed in Sec. V, the pion source-sink separation $T$ is systematically increased, while holding the time $t=\frac{T}{2}$ fixed at which the currents are inserted. Figure 4 shows the calculation of the operator $\bra{\pi(p)}\mathcal{J}_{i}(x_{0}+\xi)\mathcal{J}_{j}(x_{0})\ket{\pi(p)}$, where $\mathcal{J}_{i}$ and $\mathcal{J}_{j}$ are the antisymmetric V-A currents discussed in Sec. IV.

We perform a correlated fit to the jackknife ensemble ratios of 4pt to 2pt functions for a given momentum $p$ and spatial separation $\xi$ between the currents. In order to extract the desired matrix element from fits to our data, we assume the following single exponential form for the ratios of 4pt to 2pt functions:

[TABLE]

where $\Delta_{\rm eff}$ is the effective energy gap between the ground-state and the excited states. Therefore the ratio $\frac{C_{\rm 4pt}(T)}{C_{\rm 2pt}(T)}$ will give the desired matrix element $\bra{\pi(p)}\mathcal{J}_{i}(x_{0}+\xi)\mathcal{J}_{j}(x_{0})\ket{\pi(p)},$ up to an additional amplitude obtained from the fit to the 2pt function in the asymptotic limit of large $T$. Given the symmetries engineered into our calculation, namely the current always being inserted at $t=\frac{T}{2}$ and the source/sink interpolators being created in an identical manner, one can further assume that the excited state contamination is equal on the source and sink sides of the current.

In Figs. (3) and  (3) we present representative fits according to Eq. (34) applied to the jackknife ensemble ratio of 4pt ($C_{4{\rm pt}}$) to 2pt ($C_{2{\rm pt}}$) correlation functions as a function of source-sink separation $T$, for $\vec{\xi}=\left[0,0,a\right]$ and momenta along the $z$-direction $p=0.610$ GeV and $p=1.525$ GeV, respectively. For $p=0.610$ GeV, the largest source-sink separation we use is $T=22a$. As the signal-to-noise ratio is significantly reduced for states with large momentum, we attempt to limit additional noise in our extracted matrix elements by considering smaller separations between the currents for states with large momenta. For instance, we limit the largest source-sink separation to be $T=16a$ for $p=1.525$ GeV. As seen from Fig. (3), the ratio of the correlation function for $p=0.610$ GeV has reasonable signal for almost all the source-sink separations. However, as expected for $p=1.525$ GeV, after $T=12a$, the lattice matrix elements become very noisy and have no effect on the fit. Reasonable statistical signal in such a relatively small window of source-sink separations, compared to lower momenta data, might be a cause for concern. However it is worth remembering that with such a coarse lattice spacing ($a=0.127$ fm), $T=12a$ is sufficiently large ($\sim 1.525$ fm) to minimize any excited-state contamination. We present values of the fit parameters in Table 2.

With the fit to the data for momenta along the $z$-direction in the range $p\in\{0.610-1.525\}$ GeV and current separations $\lvert\xi\rvert\leq 4a$ in the $z$-direction, we obtain the matrix elements shown in Fig. 5. As discussed in Sec. IV, we only include $\lvert\xi\rvert\leq 4a$ in our analysis so that $\xi$ is sufficiently smaller than $\frac{1}{\Lambda_{\rm QCD}}$ and thereby ensuring the factorization of Eq. (LABEL:eq:factor).

For the lowest-order kernel, we use the following simple relation derived in Sec. IV

[TABLE]

to extract pion valence distribution $q^{\pi}_{\rm v}(x)$, by fitting the antisymmetric V-A current matrix elements $\sigma(p\cdot\xi,\xi^{2})$.

For a proper extraction of the PDF and comparison to global fit results, one would need to extend the LO matching formula in Eq. (VI) to include a higher order matching kernel between LCSs and PDFs. A next-to-leading-order (NLO) matching kernel would include $O(\alpha_{s})$ logarithmic $\xi^{2}$ and constant corrections. The logarithmic terms contain the scale dependent DGLAP evolution. The constant terms contain the information on renormalization of the lattice QCD matrix element and the partonic PDFs which leads to the scheme dependence. There also can exist higher twist effects which contaminate the results without sufficiently small $\xi^{2}$. Finally there exist potential discretization errors from the small separation size, as well as rotational symmetry breaking effects as observed in Bali:2018spj .

The extraction of the PDF using Eq. (VI) from lattice calculated data constitutes an ill-posed inverse problem. Lattice data will always be discretized and in a limited range of $\omega$. As demonstrated in Karpie:2019eiq , a naïve discretized inverse cosine transform would introduce numerical artifacts into the PDF. Solutions to this inverse problem require additional information or constraints. In the global fitting community, additional information is given in the form of smooth physically motivated functional forms as described below. PDFs extracted using this technique have been successfully shown to describe different physical processes, thereby assuring the universality of the nonperturbative PDFs. Of importance, it is known that the valence distributions of nucleon and pion are smooth functions of $x$ in the region $0<x<1$. In the spirit of the functional forms used in global fits of PDFs, we insert

[TABLE]

into Eq. (VI) and numerically perform the integration, where $N$ is the normalization such that

[TABLE]

With the limitations related to $\xi^{2}$ corrections in mind, in this preliminary calculation, we use the numerical fitting program ROOT ROOT to fit bootstrap samples of the V-A matrix elements and obtain a LO $q^{\pi}_{\rm v}(x)$-distribution. The uncertainty band in the fit has been obtained from the fit results of the bootstrap samples. With the various sources of $\xi^{2}$ corrections not taken into account, we did not expect the matrix elements as a function of Ioffe time to fall upon a single curve - consequently the $\chi^{2}/{\rm d.o.f.}$ was close to 2.2.

We find that the term $\rho\sqrt{x}$ in Eq. (36) has no effect in the fit, as $\rho\simeq 0$. A similar zero-value for $\rho$ was also found in Ref. Chen:2016sno and other global fits to experimental data. We therefore adopt a simpler functional form for the PDF in our calculation

[TABLE]

where the beta functions in the denominator ensure the normalization condition in Eq. (37) is met. In Eq. (38), the $(1-x)^{\beta}$ allows a smooth interpolation to zero as $x\to 1$ and is inspired by the counting rule of perturbative QCD. The $x^{\alpha}$ term is motivated by the behavior predicted by Regge theory at small $x$. One could interpolate these two limits using a polynomial of $x$. However, due to present statistics and small range of $\xi$, we cannot quantitatively distinguish between different choices of polynomials. Therefore, we use the widely adopted phenomenologically motivated functional form of pion valence PDF in Eq. (36). We set the following physically motivated and relaxed constraints

[TABLE]

The fit to the lattice QCD data using the LO kernel in Eq. (VI) and the functional form of PDF in Eq. (38) is shown in Fig. 5 with the fit parameters,

[TABLE]

The extracted PDF from this fit is shown in Fig. (6) where the values of the fit parameters are indicated. We also show the $xq^{\pi}_{\rm v}(x)$-distribution in Fig. (6). The perturbative kernel fixes the value of the integral in Eq. (VI) to be $\frac{1}{\pi^{2}}$ at $\omega=0$ for any value of $x$, therefore the fitted value of $T_{1}\left(\omega,\xi^{2}\right)$ has zero uncertainty at this point.

VII Comparison with Other Determinations

This first exploratory lattice QCD calculation of the pion PDF using spatially-separated current-current correlation function is performed at a relatively heavy pion mass ($m_{\pi}\simeq 416\text{ MeV}$). This calculation must be repeated on several other lattice ensembles to determine the pion mass dependence, quantify lattice artifacts such as finite lattice spacing and finite volume Briceno:2018lfj corrections and obtain the PDF in the continuum limit. As mentioned earlier, extending the perturbative calculation beyond LO will not only lead to a more reliable extraction of the PDF, but also an understanding of power corrections and higher twist effects. A NLO matching kernel will give control over the corrections in $\xi$, both from DGLAP and higher twist effects. This calculation was performed on a fairly coarse lattice with a large minimum $\xi$, and in the future these corrections will need to be taken into account. While such calculations are underway and will be presented in a future work, the limitations in our current extraction of the pion valence PDF do not preclude comparison with global fits, two different model calculations and recent lattice calculations of pion valence quasi-distribution.

For a comparison with the LO extraction of $q^{\pi}_{\rm v}(x)$ from Drell-Yan experimental data in Ref. Conway:1989fs , we evolve our lattice QCD determination of the PDF in LO to an evolution scale of $\mu^{2}=27$ GeV2 starting from initial scale of $\mu_{0}^{2}=1$ GeV2. With only a LO matching kernel, the initial scale $\mu_{0}$ is chosen to be comparable to the $\frac{1}{\xi}$’s used in this calculation, but not low enough for perturbation theory to be doubted. With a NLO matching kernel, there will exist an explicit relationship between the scales $\xi$ and $\mu_{0}$ from the logarithmic terms. After the evolution, a shift in the peak of the $xq^{\pi}_{\rm v}(x)$-distribution toward smaller values of $x$ and a more convex-up behavior of the distribution near $x=1$ is seen as expected in our calculation. From the fit parameters in Eq. (38) ($\alpha=-0.34(31)$, $\beta=1.93(68)$, and $\gamma=3.05(2.50)$ at the initial scale), it is seen that this lattice QCD calculation of $q^{\pi}_{\rm v}(x)$ is in agreement within uncertainty with the analysis in Ref. Aicher:2010cb , where the authors included next-to-leading-logarithmic threshold soft-gluon resummation effects in the calculation of the Drell-Yan cross section. The large-$x$ behavior is statistically consistent with the expectation based on perturbative QCD Farrar:1979aw ; Berger:1979du ; Brodsky:1994kg but of course with large uncertainty. In contrast, the large-$x$ behavior of this calculation has about $\sim\!\!1\sigma$ difference from the two other NLO fits Wijesooriya:2005ir ; Barry:2018ort which obtained a harder $(1-x)$ fall-off of the pion valence distribution.

It is seen in Fig. (7) that the large-$x$ behavior of this calculation is statistically consistent with the Dyson-Schwinger model prediction Chen:2016sno labeled as “DSE” in the momentum fraction region $x>0.7$. On the other hand, this lattice QCD calculation of $q^{\pi}_{\rm v}(x)$ is in statistical agreement with the light-front holographic QCD model calculation labeled as “LFHQCD” in the region $x<0.5$, but shows a slightly softer fall-off at large-$x$ in its central value. As mentioned earlier, in a future calculation, when all the systematics of this lattice QCD calculation are to be well understood and controlled in a proper way, the first-principles determination of large-$x$ behavior of pion PDF such as this one can shed light for understanding different approximations used in an array of model calculations.

Even with the limitations mentioned above, this first exploratory lattice QCD determination of $q^{\pi}_{\rm v}(x)$ using LCSs provides encouraging results and shows that this method has the potential to capture the essential dynamics dictating the behavior of hadron PDFs. Of notable interest, our calculation of $xq^{\pi}_{\rm v}(x)$ , shown in Figs. 6 and 7, illustrates a peak of the distribution in a region $x<0.50$ at any scale $\mu^{2}$. This is consistent with all the global analyses of the pion valence distribution, wherein $xq^{\pi}_{\rm v}(x)$ is peaked below $x=0.50$. The readers are referred to Chen:2018fwa and Karthik:2018wmj for recent other lattice calculations of the pion valence quasi-distribution.

VIII Summary and Outlook

We have presented the first lattice QCD calculation of the pion valence distribution using a spatially-separated vector and axial-vector current combination. We have emphasized that the spatial separation $\xi$ between the currents is a well-defined quantity in the good lattice cross sections method and plays a role analogous to capturing the correct collision dynamics in a hard scattering process and ensures the validity of factorization to obtain parton distribution functions. In this exploratory calculation, we have considered a leading-order perturbative kernel to obtain the nonperturbative valence PDF of the pion though factorization of the good lattice cross section of this vector and axial-vector matrix element. A similar calculation on other lattice ensembles is in progress to determine the pion mass dependence, quantify lattice artifacts and obtain the PDF in the continuum limit. Such a calculation of the lattice QCD matrix elements in the continuum limit, and therefore a more reliable extraction of PDF will be presented in future work with the next-to-leading-order perturbative matching kernel incorporated to understand the corrections in $\xi$ and higher twist effects. Moreover, using this most general approach, other good lattice cross sections with different current combinations will give information on different types of PDFs. Within the limitations of the present calculation, however, we would like to emphasize that the good statistical agreement between the PDF extracted here, through only leading-order factorization and the fits to the experimental data, is very encouraging and shows that this method has the potential to complement the well-established and modern state-of-art global fits of PDFs. Upon further investigation and refinement of our methodology, our lattice QCD results can be a subset of those used in future global analyses.

Acknowledgments

R.S.S., J.K., and C.E. give special thanks to Robert G. Edwards whose help has been tremendous to set up the numerical calculation. We thank Raúl A. Briceño, Martha Constantinou, Nikhil Karthik, Luka Leskovec, Tianbo Liu, Ying-Qing Ma, Nobuo Sato, Yi-Bo Yang, Jian-Hui Zhang who provided insights and expertise that greatly assisted this research.

This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Contract No. DE-AC05-06OR23177, within the framework of the TMD Collaboration. We acknowledge the facilities of the USQCD Collaboration used for this research in part, which are funded by the Office of Science of the U.S. Department of Energy. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. J.K. and C.E. are supported in part by the U.S. Department of Energy under contract DE-FG02-04ER41302 and Department of Energy Office of Science Graduate Student Research fellowships, through the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. K.O. acknowledges support in part by the U.S. Department of Energy through Grant Number DE- FG02-04ER41302, by STFC consolidated grant ST/P000681/1, and the hospitality from DAMTP and Clare Hall at Cambridge University.