A Concise Review of Gradient Models in Mechanics and Physics

The various mathematical models developed over the years to interpret the behavior of materials and corresponding processes they undergo were based on observations and experiments made at that time. Classical laws for solids (Hooke) and fluids (Navier–Stokes) form the basis of current technology. The discovery of new phenomena with the aid of newly developed experimental probes have led to various modifications of these laws, especially at small scales. The emergence of nanotechnology is ultimately connected with the design of novel tools for observation and measurements, as well as with the development of new methods and approaches for quantification and understanding. This paper first reviews the author's previously developed weakly non-local or gradient models for elasticity, diffusion, and plasticity. It then proposes a similar extension for fluids and electrodynamics. Finally, it suggests a gradient modification of Newton's law of gravity, with a possible connection to the strong force of elementary particle physics.


INTRODUCTION
In a recent chapter in Advances of Applied Mechanics [1], an extensive review of the author's internal length gradient (ILG) mechanics framework is presented. It is based on the assignment of internal lengths/ILs (associated with the local geometry/topology of material nano/microstructures) as scalar multipliers of extra Laplacian terms that are introduced into classical constitutive equations to account for heterogeneity effects and weak non-locality. Related background work for this framework can be found in the references quoted therein, as well as in earlier published articles by the author and his coworkers on gradient elasticity/plasticity [2,3] and higher order diffusion/heat conduction [4,5].
The motivation for the ILG framework goes back to van der Waals work on liquid/vapor interfaces, as extended by Landau's school through the introduction of a gradient-dependent order parameter and, more recently, adopted in Cahn-Hilliard theory of spinodal decomposition [6,7]. In this connection, it is also pointed out that the plethora of more recent phenomenological phase field models, as well as subatomic generalized functional density theories, are essentially based on the introduction of higher order gradients in the free energy of the system under consideration.
At very small scales, mechanical and chemical effects are often equipresent, and an extended chemomechanical ILG framework is necessary to consider higher order IL couplings, as suggested in Aifantis [1]. In view of the fact that mechanical and chemical ILs are introduced as scalar multipliers of corresponding Laplacian terms, it turns out that such coupled chemomechanical formulation is appealing and robust. Since in models of mathematical biology [8,9], cells are commonly represented by scalar concentration fields (i.e., in the same way as chemical species), the formulation could be easily adapted for the description of higher order couplings between mechanical and biochemical ILs. Such an extended ILG mechanics framework, including synergistic effects between mechanical and chemical or biological ILs, can be employed to consider chemomechanical instabilities in lithiumion battery/LiB anodes [10] and biomechanical instabilities in brain tumors [11]. In this connection, it should be pointed out that the Laplacian of the elastic Hookean stress in the author's model of gradient elasticity/GradEla [3,5] enters in the same way as in Murray's treatment of morphogenesis in biology [8,9]. Similarly, the Laplacians of the immotile/motile cancer cell densities in the go or grow/GoG model [11] were introduced in the same way as in the Walgraef-Aifantis/W-A model of immobile/mobile dislocation densities for persistent slip bands [12].

ILG IN SOLIDS: ELASTICITY/PLASTICITY/DIFFUSION AND CHEMOELASTICITY
For elastic deformations, the term ℓ 2 ε ∇ 2 [λε mm δ ij +2Gε ij ]-where ℓ ε denotes elastic IL, ε ij is the elastic strain, and (λ, G) are the Lamé constants-is incorporated into classical Hooke's law σ ij = λε mm δ ij + 2Gε ij , which now reads [1,3] Similarly, the term ℓ 2 p ∇ 2 γ pwhere ℓ p denotes plastic IL and γ p = γ p dt;γ p = 2ε p ijε p ij is the equivalent plastic strain with ε p ij denoting the plastic strain tensor and a superimposed "·" denotes time derivative-is introduced in the classical von Mises yield condition . For diffusion problems, the IL enters through the additional term ℓ 2 d ∇ 2 j i -where ℓ d is a diffusional internal length and j i denotes the diffusion flux-which generalizes the classical Fick's lawρ = D∇ 2 ρ by the expressionρ = D∇ 2 ρ − c∇ 4 ρ; c = ℓ 2 d D, in a manner formally similar to the Cahn-Hilliard theory [7] of spinodal decomposition. For collective dislocation phenomena, the IL enters through the extra Laplacian terms D i,m ∇ 2 ρ i,m , where ρ i,m denote (immobile, mobile) dislocation densities and D i,m "effective" diffusion-like transport coefficients [12], which are absent from conventional dislocation dynamics models.
For coupled elasto-diffusion processes, higher order IL couplings have not been adequately considered with the exception of some recent works on chemomechanical damage of LiB anodes and propagation of lithiation fronts. Colossal mechanical stresses develop during lithiation, and the following gradient enhanced stress-assisted diffusion equation may be used to model Li intercalation In the above expression, ρ is the Li concentration, σ 0 h denotes the hydrostatic component of mechanical stress, and (D, N, M) are phenomenological coefficients, whereas ℓ ρ denotes a diffusional internal length entering in an analogous (slightly different) way as ℓ d above. It is noted that Equation (1) is an IL extension of the stress-assisted diffusion equation earlier used for hydrogen embrittlement [13]. It is an uncoupled equation, since σ 0 h is assumed to be determined by the solution of a separate elasticity problem.
Higher order elasto-diffusion couplings were considered more recently in Tsagrakis and Aifantis [14] to model the propagation of lithiation fronts. The corresponding constitutive equations for the stress (σ ) and the chemical potential (µ) are given by Here again, ρ is a normalized Li concentration, ε is the strain, (λ, G) denote the Lamé constants, and (M, ) are chemomechanical phenomenological coefficients. The quantities µ 0 (reference value of µ) and α denote chemical parameters of ideal solution theory, whereas ℓ ε and κ denote the mechanical and diffusional internal lengths, respectively.

ILG IN RHEOLOGY: NEWTONIAN AND COMPLEX FLUIDS
In this section, we suggest various possibilities for a gradient enhancement of constitutive equations used in the fluid mechanics and rheology communities. In this connection, it is pointed out that (following the author's work on gradient theory), a number of such generalizations have already been proposed in these communities [15][16][17][18][19]. In the spirit of the ILG formulation, such type of generalizations can readily be deduced by formally replacing the local fields for the fluid density ρ, stretching tensor respectively. Another possibility is to include the Laplacian of the viscoelastic stress , as proposed in the diffusive Johnson-Segalman (DJS) model, i.e., with p denoting pressure, µ being the solvent's shear viscosity [not to be confused with the same symbol used for the chemical potential in Equation (2)], and the second equation standing for quasi-static equilibrium. The newly introduced viscoelastic stress obeys the following evolution equation where τ denotes relaxation time, µ * is the micelle (polymer-like) viscosity, D is a diffusion-like coefficient, and a superimposed " o " denotes corotational time derivative.
On returning to the topic of an appropriate generalization of the Navier-Stokes (N-S) equations for incompressible fluids, i.e., of the constitutive equation T = −p1 + 2µD, we can propose the following gradient model where ℓ T and ℓ D denote internal lengths associated with stress and strain rate inhomogeneities. On assuming that ℓ T can be neglected and introducing Equation (5) in the equation of momentum balance ρv = divT, we obtain the following gradient generalization of the N-S equations where = ∇ 2 and 2 = ∇ 4 denote the Laplacian and biharmonic operators, respectively [16]. A slightly generalized model to consider turbulence, reads where the α parameter [not to be connected with the same symbol used in Equation (2)] denotes a statistical correlation length and ∇ D = D +DW−WD denotes the usual Jaumann rate. Steady-state solutions of Equation (6) may be determined by employing the operator split method (or the use of Ru-Aifantis theorem [20]) utilized to eliminate singularities from dislocation lines and crack tips in the theory of gradient elasticity (see also [1]). This same procedure leads to the cancelation of singularities in typical fluid flow calculations involving immersed objects. It turns out, for example, that the resulting gradient Oseen tensor where r i denotes the position vector and r its magnitude, reads This expression is essentially the same with the exponential regularization of the Green's tensor emerged in gradient elasticity/GradEla and used to rederive non-singular expressions for the stresses and strains in dislocation lines and crack tips.
More details on such gradient fluids aspects can be found in Giusteri and Fried [21], the authors of which it seems not to have been aware of analogous developments in gradient elasticity.

ILG IN OTHER DISCIPLINES AND SCALES
In this section, we summarize the applicability of the ILG framework to other disciplines and scales ranging from earth scales to quantum scales.

ILG in Geology
Some initial work on introducing internal lengths and Laplacians of strain in geomodels has been published by the author and coworkers to model shear banding and related instability phenomena in soils/rocks and snow/ice (see, for example [22][23][24][25][26][27][28][29][30][31][32][33]). Various types of gradient-dependent constitutive equations for such classes of geomaterials have also been proposed and elaborated upon in detail by many other authors. This was due to the fact that the Laplacian regularizes unstable behavior in the geomaterial's softening regime and allows for the determination of shear band thickness/spacing, eliminates meshsize dependence, and establishes convergence of corresponding finite element calculations [2]. The popularization of the approach in the geomechanics community is mainly due to the follow-up works by Vardoulakis and collaborators for soils, as well as de Borst and collaborators for concrete. These are too many to mention here, and most of them can be found in the web.
In connection with the above, it is worth noting that the Walgraef-Aifantis/W-A model for dislocation patterning has recently been used by Ord and Hobbs [34] to interpret fracture patterns in frictional, cohesive, granular materials. Their article was one contribution of 17 to a Theme Issue "Patterns in our planet: Applications of multi-scale non-equilibrium thermodynamics to Earth-system science".

ILG in Electrodynamics
The inclusion of higher order gradients in deforming materials under the action of electromagnetic fields has also become very popular in recent years due to emerging applications and novel design of piezoelectric (induction of electricity due to applied stress/strain) and flexoelectric (induction of electricity due to strain gradients) components. The number of written articles is prohibitive for mentioning them here, and we only refer to recent publications by the author and coworkers [35][36][37], as well as to the bibliography listed there for recent literature related to size effects in micro-/nanoelectromechanical systems.
In relation to the issue of eliminating singularities and introducing screening effects (e.g., Debye screening) in the electric field, the following gradient modification of Coulomb's law of electrostatics has been used [see, for example, [38] where, in addition, a fractional generalization of Debye screening is also discussed], as follows: In this equation, is the electrostatic potential [E(r) = −∇ (r) ; E(r) is the electric field], ρ(r) denotes the charge density, ε 0 is the vacuum permittivity, r D is the Debye screening distance, and r denotes the spatial radial coordinate. The classical Coulomb's potential for a point charge of strength Q has the form (r) = Q/4πε 0 r, while its Debye-screened counterpart obtained from Equation (10) In concluding this discussion on gradient electrodynamics, an interesting possibility is outlined below for a gradient extension of the MacCullagh-Fitzgerald formal analogy between elastic deformation and light transmission/reflection/refraction in transparent media (also discussed briefly by Sommerfeld in his lectures), as pointed out in Truesdell and Toupin [39]. On assuming first that the ether behaves as an elastic medium with its stress T depending linearly on rotations ω(instead of strains), we can write where u denotes displacement, ρ denotes density, and k is a constant. These lead to the equation k curl curl u + ρü = 0 and by formally setting k curl ü → α E and ρu → αβ (with E/B denoting electric/magnetic fields and α being an arbitrary proportionality constant), we arrive at Maxwell's equations [39] ∂B ∂t where the identities divcurlu = 0 and curl (∂u/∂t) − ∂(curl u)/∂t = 0 were used, along with the following identification of the various coefficients (β = kε 0 , µ 0 = ρ/β, with β being an arbitrary constant). The symbols (ε 0 , µ 0 ) are the fundamental vacuum constants of electromagnetism with ε 0 µ 0 = 1/c 2 and c denoting the speed of light in vacuum. By repeating the same procedure as above but replacing Equation (12) 1 for the ether's elastic stress with its gradient counterpart we arrive at the following generalization of Maxwell's equations It is noted that for electrostatics (with the usual assumption that the electric field E is proportional to a potential gradient ∇ ), Podolsky's non-quantum equation for the electric potential,

ILG in Atomistics and Quantum Mechanics
We conclude this section on the applicability of the ILG framework to various disciplines and scales by focusing on two specific topics: a possible gradient generalization of the molecular dynamics stress, and an analogous generalization of the quantum mechanical stress. In this connection, it is noted that the following expressions were proposed for these stresses: 1 Podolsky has derived a generalization of Maxwell's equations through a variational principle, leading also to the appearance of ∇ 2 B in addition to∇ 2 E. This is also possible through the aforementioned analogy by replacing u with u−ℓ 2 ∇ 2 u.
in Zimmerman et al. [42], and (17) in Maranganti and Sharma [43], with the various symbols having their usual meaning [42,43]. The formal similarity between these two expressions and their resemblance to the virial stress and other statistical stress measures is striking. However, the problem to connect such discrete "microscopic" stress measures with the continuum "macroscopic" measure of Cauchy stress in a "seamless" way is a challenging issue. A gradient generalization of the force fields f ij in Equation (16) and the interaction potential U ij in Equation (17) may turn out to be quite useful in this respect. Among other things, it could also naturally introduce screening distances and eliminate associated singularities.
The effect of strain ε on the electronic structure has been described [44] through the equations where ψ(r) denotes the wavefunction, [C] is the Hookean (generally anisotropic) elasticity matrix, a c is the so-called deformation potential constant, and the rest of the symbols have their usual quantum mechanical meaning [44,45]. This is an uncoupled framework where strain can affect the electronic state but not vice versa. A generalization to account for the inverse effect on strain due to changes in the quantum field through the wavefunction ψ(r) has already been proposed as follows [45]: (E c − ℏ 2 2m * ∇ 2 )ψ(r) + a c tr(ε)ψ(r) = Eψ(r) and (19) where K is the isotropic bulk elastic modulus. A possible gradient modification is then to replace ε with its gradient counterpart ε − ℓ 2 ε ∇ 2 ε, and this formal generalization may be of interest to further explore.

ILG MODIFICATION OF NEWTON'S GRAVITATION LAW
In this final section, we venture a gradient generalization of Newton's Law, which allows for the corresponding gravitational force to attain values larger than the electromagnetic force and even reach the levels of the nuclear/strong force which keeps matter together 2 .
We begin with the following nonlocal generalization of the gravitational force f in its component form (f i ) where G ij (r − r ′ ) is a non-local interaction kernel, and F j is the classical Newton's force. By Fourier transforming Equation [26], Taylor series expanding up to the second-order term and inverting, the following differential equation is obtained for the above listed expression of the generalized force f where k = |k| denotes wave vector,G ij is the Fourier transform of G ij , and ℓ is an internal length, with δ ij appearing due to the assumed isotropy/spherical symmetry [In general the sign in front of the Laplacian term in Equation (21) may be positive or negative depending on the sign of d 2G ij (0)/dk 2 of the second order term in the Taylor expansion without assuming its absolute due to stability considerations]. Such a formal derivation can also be established by viewing the two point masses M 0 and M in the classical Newton's Law, as being distributed and bounded by spheres of finite radii. By considering, for example, the reference mass M 0 (M 0 = i m i ) as being distributed within a sphere of radius R 0 , summing up the interactions of each point mass m i (located at distance r i from the center of the sphere where r i = 0) with the other distant mass M (viewed as a point mass), and expanding in Taylor series the density ρ(r i ) around ρ 0 = ρ(0) by keeping terms up to the second order, we obtain the following relationship where e R = R/R denotes the unit vector along the line connecting the center of the reference mass M 0 with the second distant point mass M. On setting V ρ 0 dV = M 0 , we then have On assuming a radial dependence of f and F[f = f (r) e r ; F = F e r , F = A/r 2 , with A denoting Newton's classical gravitational proportionality parameter and setting e R = e r ], we can readily solve the scalar radial counterpart of Equation (21) in spherical coordinates, i.e.
by also requiring that F → 0 as r → ∞. The result is where B is a new parameter to evaluate in connection with experiments. It is noted that the above expression of Equation [30] reduces to Newton's classical force F N = A/r 2 as r → ∞ and to the expression F SF = AB/r 2 as r → 0. By adjusting the value of the new parameter B (B >> 1), we can attain values of the nuclear and strong force for F SF , as outlined below. First, we note that the internal length parameter (ℓ) can be identified with de Broglie relativistic length, the Compton length, the Planck length, or the Schwarzschild distance, according to the configuration at hand. If, for example, ℓ is identified with de Broglie relativistic length used in the Vayenas [46] rotating neutrino model (RNM), it follows that ℓ =h/γ m 0 c = 6.31 × 10 −16 m, whereh denotes the Planck constant, c is the speed of light, and m 0 denotes the rest mass of neutrino, whereas G is the classical Newton's gravitational constant (not to be confused with the same symbol earlier used for the shear modulus), and γ is the Lorentz factor (γ = 1/ 1 − (v/c) 2 ; with v denoting particle speed).
To proceed further, we adopt next the Vayenas and coworkers RNM method [46,47] as applied to the nucleus (proton or neutron), but we use the gradient enhanced Newton's Law as given by Equation (24) with B = 0 instead of the classical Newton's Law expression (B = 0). To this end, we shall equate f in Equation (24) with the centrifugal force F C = γ m 0 c 2 /r. Here, r denotes the radius of the nucleus modeled by the three rotating neutrinos whose total relativistic mass is m N = 3γ m 0 . For a proton nucleus (m N ≡ m p ), an estimate of γ can be obtained by equating the proton energy m p c 2 with the one corresponding to the relativistic neutrino mass. This gives the value of γ = m p 3m 0 , which, according to experimental measurements for m p (= 9.38 × 10 8 eV/c 2 = 1.67 × 10 −27 kg) and m 0 (= 0.04 eV/c 2 = 7.13 × 10 −38 kg) turns out to be equal to 7.818 × 10 9 .
Having such a value of γ available, we can make effective use of the equality between gravitational and centrifugal forces, as proposed in the RNM approach, to deduce the relationship where the factor √ 3 rises by considering the resultant gravitational force F R = f / √ 3. One possibility for the constant A is to assume that it is given by the expression A = Gm 2 0 γ 2 to account for relativistic effects (during the interaction of each pair of the rotating neutrinos in the assumed RNM configuration) in the same manner as in the expression for the centrifugal force F C (where the square is now introduced due to the fact that two neutrinos are involved). The above relationship (with ℓ identified with de Broglie's relativistic length, as discussed in the previous paragraph, i.e., ℓ = 6.31 × 10 −16 m) gives The resulting value of F R is i.e., the value of the strong force obtained for the RNM configuration [46,47] using an entirely different approach. In that approach, Equation (25) with B = 0 was used with A = Gm 2 0 γ 6 , giving a value for γ = 3 1/12 m 1/3 pl m −1/3 0 = 7.167 × 10 9 , where m pl is the Planck mass (m pl = √ ℏc/G), and the value of m 0 was taken as m 0 = 0.0436 eV/c 2 . Since γ = m p /3m 0 , this gives m p = 9.38 × 10 8 eV/c 2 , i.e. the same value as the one used in the previous paragraph by properly adjusting the parameters (A, B), as well as by identifying the internal length parameter ℓ with de Broglie relativistic length. Other choices of (A, B, ℓ ) are possible not only for the RNM configuration at hand but also other more complex geometric models for elementary particles represented by several neutrinos where a potential is convenient to use. This topic will be pursued independently in a future article.

CONCLUSIONS
A concise review of gradient models (across scales, materials, and processes) was provided based on the author's ILG approach. As a result, earlier references on generalized continuum mechanics and recent contributions on gradient and non-local theories were not discussed due to space limitation. For solids, one should single out the contributions of Eringen [48], Fleck and Hutchinson [49,50], Gurtin and Anand [51], Gao et al. [52], Nix and Gao [53], de Borst et al. [54,55], Geers et al. [56], Peerlings et al. [57], Willis et al. [58], Aifantis and Willis [59], and Polizzotto [60,61]. Many more are included in a most recent and detailed article by Voyiadjis and Song [62] focusing on gradient plasticity. Gradients in fluid and granular flows were considered most recently by Goddard [63,64]. For additional recent developments on granular flow, one may also consult references [65][66][67][68][69], while for internal length interpretations based on kinetic theory, one may consult [70]. After writing this article, it came to the attention of the author that an expression similar to that derived herein and given by Equation (24) was also proposed on rather intuitive grounds by Fischbach et al. [71] in an effort to reinterpret existing measurements on earth's gravity (see also [72]). The values of their constants were entirely different than ours, as they used it for a reanalysis of the Eötvös experiment on Earth's gravitational field. There has been a vast literature on this expression, subsequently referred to as the "fifth force, " which we will discuss in a forthcoming publication, as this is beyond the scope of the present review.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

ACKNOWLEDGMENTS
Useful discussions with my daughter Katerina Aifantis and my student Kostas Parisis along with my former university classmate Constantinos Vayenas, as well as the invitation of Joshua Dijksman to contribute to this volume are gratefully acknowledged. The work was greatly benefited from interactions within the RISE/FRAMED project no. 734485 (https://cordis. europa.eu/project/rcn/207050/factsheet/en) for which Aristotle University of Thessaloniki (AUTh) acts as coordinator. In this connection, thanks are extended to all beneficiary nodes and international partners of FRAMED. The gradient fluids section is a topic of the RISE/ATM2BT project no. 824022 (https:// cordis.europa.eu/project/rcn/219192/factsheet/en) for which EU support is also acknowledged. This section was included in anticipation of a follow-up joint work between AUTh, Akita Univeristy/Japan, and Aston University/UK (which acts as project's ATM2BT coordinator). More details on this minireview can be found in a forthcoming arXiv preprint.