Revisiting Lorentz violation in Hořava gravity

# Revisiting Lorentz violation in Hořava gravity

Andrew Coates Theoretical Astrophysics, IAAT, University of Tübingen, Geschwister-Scholl-Platz, 72074, Tübingen, Germany School of Mathematical Sciences, University of Nottingham, University Park, Nottingham, NG7 2RD, United Kingdom    Charles Melby-Thompson Institut für Theoretische Physik und Astrophysik, Julius-Maximilians-Universität Würzburg, Am Hubland, 97074 Würzburg, Germany Department of Physics, Fudan University, 220 Handan Road, 200433 Shanghai, China    Shinji Mukohyama Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
July 13, 2019
###### Abstract

In the context of Hořava gravity, the most promising known scenarios to recover Lorentz invariance at low energy are the possibilities that (1) the renormalization group flow of the system leads to emergent infrared Lorentz invariance, and (2) that supersymmetry protects infrared Lorentz invariance. A third scenario proposes that a classically Lorentz invariant matter sector with controlled quantum corrections may simply co-exist with Hořava gravity under certain conditions. However, for non-projectable Hořava gravity in dimensions it is known that, in the absence of additional structures, this mechanism is spoiled by unexpected power-law divergences. We confirm this same result in the projectable version of the theory by employing the recently found gauge-fixing term that renders the shift and graviton propagators regular. We show that the problem persists for all dimensions , and that the degree of fine tuning in squared sound speeds between a gauge field and a scalar field increases with . In particular, this difference in the zero external momentum limit is proportional to for , where is the ultraviolet momentum cutoff for loop integrals, while the power-law divergences are absent for and . These results suggest that not only the gravity sector but also the matter sector should exhibit a transition to Lifshitz scaling above some scale, and that there should not be a large separation between the transition scales in the gravity and matter sectors. We close with a discussion of other more promising scenarios, including emergent Lorentz invariance from supersymmetry/strong dynamics, and pointing out challenges where they exist.

preprint: YITP-18-11, IPMU18-0034

## I Introduction

Despite the empirical successes of General Relativity (GR), there are reasons to believe that it is incomplete. The major stumbling block from a theoretical perspective is that quantum GR is not perturbatively renormalizable. One approach to improving its UV properties is to relinquish Lorentz invariance at short distances Horava:2009uw (), allowing the inclusion of higher order derivatives that render gravity renormalizable without violating unitarity. Renormalization group (RG) flow is organized around a fixed point with “Lifshitz scaling”, which acts anisotropically on space and time, so that higher order spatial derivatives scale the same way as second order time derivatives in the ultraviolet Horava:2009uw (); Visser:2009fg ().111The Lifshitz scaling is also the origin of interesting cosmological implications of the theory, such as a novel mechanism for generating scale-invariant perturbations Mukohyama:2009gg () and a solution to the flatness problem Bramberger:2017tid (). See Mukohyama:2010xz () for a review of cosmology based on Hořava gravity. We refer to such theories as Hořava gravity.

Lorentz violations in the matter sector, however, are very tightly constrained. For example, the relative differences among speed limits for different species of matter fields are typically constrained to be smaller than or so. Moreover, thanks to the recent multi-messenger observation of a binary neutron star merger, we now know that even the speed of gravitational waves cannot be different from that of photons by more than or so. Therefore, naturalness demands a mechanism either protecting infrared Lorentz symmetry, or causing it to emerge Kharuk:2015wga (); Afshordi:2015smm () in the infrared. If Hořava gravity or something similar is the correct description of gravity in the nature, those constraints require the RG flow of Lorentz-violating couplings to vanish rapidly in the IR Chadha:1982qq (). This may be the case if supersymmetry protects infrared Lorentz invariance GrootNibbelink:2004za (), and/or the system passes through a strongly coupled region which speeds up the otherwise logarithmic RG running of Lorentz-violating couplings toward zero Anber:2011xf (); Bednik:2013nxa (). We briefly discuss these possibilities in section IV.

There is yet another possibility that, until recently, seemed quite promising: the mechanism proposed in Pospelov:2010mp (). The basic statement of this mechanism is that, if one assumes that matter is relativistic at tree-level, then Planck suppression of interactions between gravity and matter can lead to a suppression of Lorentz-violating corrections by the ratio

 δc2c2∼(M∗Mp)2∼10−18(M∗1010GeV)2, (1)

where is the scale suppressing the higher order spatial derivatives in the gravity sector. However, it was already shown in Pospelov:2010mp () that the vector sector of the non-projectable version of Hořava gravity does not benefit from this mechanism. Instead, one finds divergences of the form

 δc2c2∼(ΛMp)2, (2)

where is the ultraviolet momentum cutoff. There, where the temporal gauge was used, this was attributed to the lack of Lifshitz scaling in the vector sector.222The authors state, however, that they verified their computations in a generalization of the gauge.

The authors of Pospelov:2010mp () proposed eliminating the power law divergences by introducing mixed derivatives into the Hořava action, in the form of terms quadratic in ( is the extrinsic curvature of the spatial slices), see also Colombo:2015yha (); Colombo:2014lta (). In the non-projectable model, however, radiative corrections in the presence of mixed derivatives is expected to give rise to terms quadratic in time derivatives of the lapse , leading to a new scalar degree of freedom that is unstable in the IR Coates:2016zvg (). The viability of the mechanism, as presented in Pospelov:2010mp () for the non-projectable theory, therefore needs a more in-depth analysis. We shall briefly discuss this point in section IV.1.

The present paper has two main goals. The first is to investigate the same mechanism in the projectable version of Hořava gravity in dimensions. The projectable theory has technical advantages over the non-projectable theory, not the least of which is that it is known to be renormalizable Barvinsky:2015kil (), while renormalizability of the non-projectable theory remains an open question. A technical but crucial difference between the two theories is that only in the projectable case are all propagators in the gravity sector regular, explicitly respect the Lifshitz scaling in the ultraviolet, and have no instantaneous modes,333In the non-projectable theory there exist modes whose momentum space correlators diverge with respect to spatial momentum even at non-zero frequency , which leads to instantaneous propagation. provided we adopt the gauge fixing of Barvinsky:2015kil (). It is therefore beneficial to revisit Lorentz violation in the matter sector due to quantum corrections in the context of the projectable version of the theory. We shall confirm that the divergence structure is robust and that the problem persists in dimensions for any .

Our second goal is to discuss possible resolutions of these problems. These fall into three classes: mixed derivative terms in the projectable Hořava gravity, supersymmetry, and strong dynamics. The mixed derivative terms proposed in Pospelov:2010mp () for the non-projectable model do not introduce a new scalar degree of freedom in the projectable theory, because there the offending term does not exist. One is instead faced with infrared instabilities, imposing strong conditions on the value of . The viability of the mechanism then reduces to the RG properties of , together with the existence of an analogue of the Vainshtein mechanism rendering the dynamics at near 1 regular. A second option — first investigated in GrootNibbelink:2004za () — is the use of supersymmetry to control Lorentz violation in the infrared, and we comment on the possibility of combining the mechanisms of GrootNibbelink:2004za (); Pospelov:2010mp (); Pospelov:2013mpa () to suppress the loop contributions of the shift variable. Finally, we briefly discuss the strong coupling mechanism studied in Anber:2011xf (); Bednik:2013nxa ().

The rest of the paper is organized as follows. Section II introduces our models, consisting of Hořava gravity coupled to a free scalar and a gauge field respectively, gives their propagators in the gauge-fixing of Barvinsky:2015kil (), and lists the integrals that appear when calculating divergent contributions to the one-loop effective action. Section III computes the one-loop corrections to the scalar and gauge limiting speeds induced by the gravitational sector. Possible solutions to the naturalness problem for Lorentz symmetry are discussed in section IV, and section V summarizes our conclusions and discusses interesting directions for future work.

## Ii Setup and notation

Our treatment of Hořava gravity follows the setup of Barvinsky:2015kil () closely. We shall only take the extrinsic curvature terms and the terms which are of order in spatial derivatives, as we are interested only in the UV contributions of gravitational degrees to loop integrals. We shall use the gauge-fixing prescription of Barvinsky:2015kil (), for which all propagators are regular. On the other hand, we needn’t worry about introducing Faddeev-Popov ghost fields, since at one loop they do not contribute to the matter sector effective action. As we consider the projectable version of the theory, we are also free to set the lapse, , to throughout the paper. We shall work in Euclidean signature, so that below has a minus sign relative to the real time action. We then add a matter sector that (before Wick rotation) is Lorentz-invariant at the classical level, and investigate Lorentz-violating quantum corrections in the zero external momentum limit.

### ii.1 Gravity action

Hořava gravity starts with a fixed foliation of spacetime by constant time slices. In terms of a coordinate system respecting the foliation, it is built from the ADM fields , , and . We work solely with the so-called “projectable” theory, where is constant on spatial slices. Its action is

 SHL=12κ2∫dtdDxN√g[KijKij−λK2+LV+Lg.f.], (3)

where , and is built from the Riemann curvature tensor of the spatial slice using terms containing spatial derivatives. The coupling constants and are dimensionless under the UV scaling.

To compute the leading divergent contributions to the matter effective action, it is enough to set the cosmological constant to zero and expand around a flat background:

 N=1,Ni=Ni,gij=δij+hij, (4)

where and are small perturbations. The gauge fixing term of Barvinsky:2015kil () has two free parameters, and , and takes the form

 Lg.f.=σFiOijFj, (5)

where

 Oij = Δ−D+2[δijΔ+ξ∂i∂j]−1, (6) Fi = ˙Ni+12σO−1ij∂khjk−λ2σO−1ij∂jh. (7)

This specific choice for the gauge fixing term decouples the shift and metric perturbations at quadratic order and renders all propagators regular. That this gauge-fixing term takes a universal form depending in a simple way on dimension can be traced to the fact that only appears in the kinetic term, which is independent of dimension and the choice of potential in the projectable theory.

### ii.2 Matter action

Following Pospelov:2010mp (), we take the matter sector to be Lorentz invariant at the HL scale, and evaluate the leading Lorentz-violating quantum corrections in the limit of small external momentum. We will consider two types of matter: a massless scalar and a gauge field, both taken to be Lorentz invariant in the UV. After Wick rotation, the scalar action takes the form

 Sϕ=12∫dD+1x[(∂tϕ−Ni∂iϕ)2+β2δij∂iϕ∂jϕ+O(hij)], (8)

while the gauge field has the action

 SA=∫dD+1x[12(Ei−NkFki)2+β24FijFij+O(hij)],Ei=DEμEiAμ,Fij=DFμFijAμ, (9)

with and . Here, we have suppressed the couplings to because, as we argue in the following sections, loops contribute only logarithmic divergences. The parameter () is the bare propagation speed of and , which we include for book-keeping purposes.

### ii.3 Propagators

Because of the particular form of the gauge fixing term, the shift propagator in Fourier space takes a simple form valid in any dimension:

 ~GijN(Pμ)=a1(δij−^pi^pj)P1(Pμ)+a2^pi^pjP2(Pμ),a1=κ2σ,a2=κ2(1+ξ)σ, (10)

where is the -momentum, , and

 P1,2(Pμ)=p2(D−1)ω2+α21,2p2D,α21=12σ,α22=(1−λ)(1+ξ)σ. (11)

For the graviton , any potential satisfying physically reasonable inequalities leads to a well-behaved Lifshitz dispersion relation, , for all modes. This leads to a Fourier-space propagator of the form

 Ghijkℓ(Pμ)=∑nM(n)ijkℓP′n(Pμ), (12)

where runs over a set of tensor structures (projecting onto transverse traceless, divergenceless vector, and two scalar modes), and the are all of the form

 P′n(Pμ)=a′nω2+α′2np2D. (13)

The gauge-fixing term is such that the mixed propagator vanishes.

For a scalar with the action (8), the Fourier space propagator is

 ~Gϕ(Pμ)=1ω2+β2p2, (14)

while a vector field in Feynman gauge with the action (9) has Fourier-space propagator

 ~GAμν(Pρ)=ηEμνω2+β2p2, (15)

where , and .

### ii.4 Loop integrals

The leading divergence of the one-loop effective action for and will be derived in the next section. They can be expressed in terms of the following loop integrals:

 J(I)1≡aI∫dD+1P(2π)D+1PI(Pμ),J(n)4≡∫dD+1P(2π)D+1P′n(Pμ),J(I)2≡aI∫dD+1P(2π)D+1ω2~Gϕ(Pμ)PI(Pμ),J(n)5≡∫dD+1P(2π)D+1ω2~Gϕ(Pμ)P′n(Pμ),J(I)3≡aI∫dD+1P(2π)D+1p2~Gϕ(Pμ)PI(Pμ),J(n)6≡∫dD+1P(2π)D+1p2~Gϕ(Pμ)P′n(Pμ),J7≡∫dD+1P(2π)D+1~Gϕ(Pμ). (16)

These integrals are regularized and evaluated in appendix A.

## Iii One-loop corrections in D spatial dimensions

Our goal is to compute the power-law divergent contributions to the shift in the limiting speeds of the fields and . This is done by computing the one-loop effective action for each field. We will see in sections III.2 and III.4 that the divergence from is logarithmic, and so the power-law divergence comes entirely from the shift . We therefore begin with loop corrections from the shift field.

### iii.1 Shift field loop corrections to the scalar action

To compute the contribution of shift loops, we may set and consider the action

 S(ϕ,Ni)=12∫dD+1x[(∂tϕ−Ni∂iϕ)2+β2δij∂iϕ∂jϕ+NiKijNj], (17)

where is the gauge-fixed kinetic term for . We are interested in the one-loop effective action for ,

 Γ1(ϕ)=12tr(logS(2)[ϕ]−logS(2)[0]), (18)

where is the fluctuation operator, with matrix components

 [S(2)[ϕ]⋅f](x)=∫dD+1yδ2S(φ)δφ(x)δφ(y)∣∣∣Ni=0f(y), (19)

stands for , and is a test function. Let denote the propagator matrix, and define the differential operator by

 G⋅S(2)[ϕ]=1+G⋅U[ϕ]. (20)

Note that . Taylor expanding in , we obtain

 Γ1(ϕ)=12tr(G⋅U[ϕ]−12G⋅U[ϕ]⋅G⋅U[ϕ]+O(ϕ3)). (21)

The components of are given by:

 ∫dd+1xdd+1yf(x)U[ϕ]ϕ(x)ϕ(y)g(y) =0, ∫dd+1xdd+1yf(x)U[ϕ]ϕ(x)Ni(y)vi(y) =∫dd+1xf(x)[∂0(vi(x)∂iϕ(x))+∂i(vi(x)∂0ϕ(x))], ∫dd+1xdd+1yvi(x)U[ϕ]Ni(x)ϕ(y)f(y) =∫dd+1xf(x)[∂0(vi(x)∂iϕ(x))+∂i(vi(x)∂0ϕ(x))], ∫dd+1xdd+1yvi(x)U[ϕ]Ni(x)Nj(y)uj(y) =∫dd+1xvi(x)uj(x)∂iϕ(x)∂jϕ(x), (22)

where , , and are test fields. Since the Green’s function vanishes, the effective action is given to quadratic order in by

 Γ(2)1(ϕ)=12tr(GijNUNiNj)−12tr(GϕUϕNiGijNUNjϕ). (23)

Consider the first trace in (23), corresponding to the diagram in figure 2. In the Lifshitz gauge fixing Barvinsky:2015kil (), the (momentum space) propagator for is a sum of two terms as shown in (10). When appears inside an integral multiplying a scalar expression, rotational invariance allows us to replace . The relevant integrals are defined in subsection II.4 and evaluated in appendix A. The final form of the divergent contribution to the effective action from the first term of (23) is therefore

 ∫dD+1xGijN(x,x)∂iϕ∂jϕ(x)=12∫dD+1xζ3δij∂iϕ(x)∂jϕ(x) (24)

with

 ζ3=δijDGijN(x,x)=2∑I=1AIJ(I)1,A1=1−1D,A2=1D. (25)

From the second term in (23) (with diagram in figure 2), we obtain a contribution to of the form

 −12∫dD+1x∫dD+1y[(−˙ϕ∂j−∂jϕ∂τ)(x)Gϕ(x,y)(˙ϕ∂i+∂iϕ∂τ+2∂i˙ϕ)(y)GijN(y,x)]. (26)

From the shift symmetry of we know that all contributions to the Lagrangian must be of the form where and are differential operators of order at least one. We are interested in the coefficients of and . (Note that by time reversal, parity, and rotational symmetry, all other contributions to the effective action are less divergent, and have more derivatives acting on ). This leaves the leading divergent contribution

 12∫dD+1x[ζ1(˙ϕ)2+ζ2(∂ϕ)2]. (27)

For we find

 ζ1=[∫dD+1y∂(x)jGϕ(x,y)∂(y)iGijN(y,x)]div=−J(2)3; (28)

note that does not contribute because its prefactor inside is transverse. For , both and contribute as

 ζ2=−[δijD∫dD+1y∂(x)τGϕ(x,y)∂(y)τGijN(y,x)]div=−2∑I=1AIJ(I)2. (29)

We thus find the divergent contribution to the effective action quadratic both in and in derivatives to be

Renormalizing so that has canonical kinetic term , we may read off the loop contribution to the infrared shift in the limiting speed:

 δc2ϕ=ζ2+ζ3+β2δZ(1)ϕ=ζ2+ζ3−β2ζ1=2∑I=1AI(J(I)1−J(I)2)+β2J(2)3=β2[(1−1D)J(1)3+(1+1D)J(2)3], (31)

where is the part of the field renormalization constant. The last equality follows from the first identity of equation (68).

### iii.2 Graviton loop corrections to scalar action

We now show that the one-loop contributions from the graviton to do not contain power-law divergences. One-loop gravitational contributions to the matter propagator come in two types, each of which is summed over the contributions of the two graviton modes. (By “graviton” we mean here the modes contributing to the correlator (12).) The first is the single-vertex loop, and the second is the two-vertex loop describing the emission and subsequent reabsorption of a virtual graviton (see figures 2 and 2).

In the first case, the loop integral involves only a factor of the graviton propagator, contributing a linear combination of integrals of the form , as defined in section II.4. These integrals are shown in Appendix A to be free from power-law divergences.

The second diagram involves in addition a factor of the scalar propagator. The vertices are a sum of and contributions, and as a result the diagrams are built from sums of the scalar and graviton propagators, weighted by the square of either the internal scalar energy or momentum. As a result, the leading divergence of these diagrams is given by a linear combination of integrals of the form and , also defined in section II.4. Again, these contributions are free from power-law divergences, as shown in Appendix A. We have thus established that the power law divergence of the one-loop correction to the infrared limiting speed of the scalar field is given by (31).

### iii.3 Shift field loop corrections to gauge field action

For shift loop contributions to the gauge field, we once again set , and consider

 S=∫dD+1x[12(Ei−NkFki)2+β24FijFij+NiKijNj]+Sgauge-fixing, (32)

where and are defined in (9), is the gauge-fixed kinetic term for , and is the same bare propagation speed as before. We adopt the Feynman gauge for , so that the propagator is given by (15), although the final result is gauge invariant.

This time we are interested in the quadratic one-loop effective action for :

 Γ(2)1(A)=12tr(G⋅U[A]−12G⋅U[A]⋅G⋅U[A]), (33)

where the operator has components

 ∫dd+1xdd+1yXμ(x)U[A]Aμ(x)Aν(y)Yν(y) =0, ∫dd+1xdd+1yXμ(x)U[A]Aμ(x)Nj(y)vj(y) =∫dd+1xXμ(x)[EkDFμFjkvj+FjkDEμEkvj+2(DEμEkFjk)vj], ∫dd+1xdd+1yvi(x)U[A]Ni(x)Aμ(y)Yμ(y) =−∫dd+1xvi(x)[EkDFμFikXμ(x)+FikDEμEkXμ(x)], ∫dd+1xdd+1yvi(x)U[A]Ni(x)Nj(y)uj(y) =∫dd+1xvi(x)uj(x)FikFjk. (34)

Here, and as in subsection II.2, and , , and are test fields. We first compute:

 (35)

where . Next, we have the contribution. Pulling this out of the analog of (26), we have

 −12tr(GN⋅UNA⋅GA⋅UAN)∣∣∣E2=12∫dD+1x∫dD+1y(EℓDFμFjℓ)(x)GAμν(x,y)(EkDFνFik)(y)GijN(y,x). (36)

As before, the dominant divergence is

 12∫dD+1xEkEℓZkℓ, (37)

where and

 ~Zkℓ =∫dD+1yDFμ(x)FjℓGAμν(x,y)DFν(y)FikGijN(y,x). (38)

By rotational invariance, , where

 ζA1=1DδkℓZkℓ =1D∫dD+1y∫dD+1yDFμ(x)FjkGAμν(x,y)DFν(y)FikGijN(y,x) =−1D∫dωdDp(2π)D+1(pjηEkν−pkηEjν)(piδνk−pkδνi)ω2+β2p2[a1(δij−^pi^pj)P1+a2^pi^pjP2] =−(1−1D)∫dωdDp(2π)D+1p2ω2+β2p2δij[a1(δij−^pi^pj)P1+a2^pi^pjP2] =−(1−1D)2∑I=1J(I)3. (39)

Similarly, the term is

 12∫dD+1x∫dD+1y(FjℓDEμEℓ)(x)GAμν(x,y)(FikDEνEk)(y)GijN(y,x), (40)

 14∫d3xFikFjℓZik,jℓ, (41)

where and

 ~Zik,jℓ=2∫dD+1yDEμ(x)EℓGAμν(x,y)DEν(y)EkGijN(y,x)∣∣∣div. (42)

Using , , and , we obtain

 ~Zik,jℓ = −2∫dωdDp(2π)D+1(ωδμℓ−δμ0pℓ)ηEμνω2+β2p2(ωδνk−pkδν0)[a1(δij−^pi^pj)P1+a2^pi^pjP2] (43) = −2∫dωdDp(2π)D+1ω2δkℓ+β2pkpℓω2+β2p2[a1(δij−^pi^pj)P1+a2^pi^pjP2].

For each tensor structure inside the integral we can replace

 pipj→1Dp2δij,pipjpkpℓ→0, (44)

where the latter comes from contracting with . Taking into account index symmetries, the final result can be written as . The contribution to from the term is

 −2∫dωdDp(2π)D+1p2(D−1)(ω2(1−1D)+β2p2(1D))(ω2+α21p2D)(ω2+β2p2)=−2[(1−1D)J(1)1−(1−2D)β2J(1)3] (45)

while the contribution from the term is

 −2D∫dωdDp(2π)D+1ω2p2(D−1)(ω2+α21p2D)(ω2+β2p2)=−2DJ(2)2. (46)

This gives the net result

 ζA2 =−2[(1−1D)J(1)1−β2(1−2D)J(1)3+1DJ(2)2]. (47)

Combining these results and the contribution from field renormalization, one obtains

 δc2A = (ζA2+ζA3)−β2ζA1=2D(J(2)1−J(2)2)+β2[2(1−2D)J(1)3+(1−1D)2∑I=1J(I)3] (48) = β2[(3−5D)J(1)3+(1+1D)J(2)3].

Once again, the last equality follows from the first identity of (68).

### iii.4 Graviton loop corrections to gauge field action

The argument given in section III.2 that the graviton modes do not contribute power law divergences to carries over to the gauge field without modification. Therefore the power law divergence of the one-loop correction to the infrared limiting speed of the gauge field is given by (48).

### iii.5 Lorentz-violating quantum corrections

To summarize, we find that the shift loop contributions (31) and (48) generate power law divergences in the difference in squared sound speeds between a gauge field and a scalar field, of the form

 δc2A−δc2ϕ=2(D−2)DJ(1)3. (49)

This vanishes only for , where vanishes because of the absence of the first term () in (10), and for . As shown in Appendix A, the integral diverges as , where is the ultraviolet momentum cutoff. Therefore, we have a naturalness problem for all . For (where ), this reduces to

 δc2A−δc2ϕ=14π2(ΛMp)2. (50)

## Iv Possible solutions

In this section we discuss some possible resolutions coming from modification of the gravity theory. Since we have shown that the naturalness problem for Lorentz-violating couplings only becomes worse for , in what follows we shall focus on .

### iv.1 Mixed derivatives

The first proposed solution to this naturalness problem, presented in the original paper, was to introduce mixed derivative terms Pospelov:2010mp (). In non-projectable Hořava gravity this means allowing terms quadratic in the following objects:

 ∇iKjk,rij≡1N(˙Rij−Nk∇kRij−Rik∇jNk−Rkj∇iNk),Ai≡12N(˙ai−Nj∇jai−aj∇iNj), (51)

where . Unfortunately, the term quadratic in the last of these,