Communications in Theoretical Physics

PAPER

A Derivation of the Entropy-Based Relativistic
Smoothed Particle Hydrodynamics by Variational
Principle
To cite this article: Philipe Mota et al 2017 Commun. Theor. Phys. 68 382

 

View the article online for updates and enhancements.

Related content
Improved kernel gradient free-smoothed
particle hydrodynamics and its applications
to heat transfer problems
Juan-Mian Lei and Xue-Ying Peng

-

Entropy-based relativistic SPH
T Kodama, C E Aguiar, T Osada et al.

-

Smoothed particle hydrodynamics
J J Monaghan

-

Recent citations
On the peripheral tube description of the
two-particle correlations in nuclear
collisions
Dan Wen et al

-

This content was downloaded from IP address 186.217.236.64 on 12/06/2019 at 15:53

https://doi.org/10.1088/0253-6102/68/3/382
http://iopscience.iop.org/article/10.1088/1674-1056/25/2/020202
http://iopscience.iop.org/article/10.1088/1674-1056/25/2/020202
http://iopscience.iop.org/article/10.1088/1674-1056/25/2/020202
http://iopscience.iop.org/article/10.1088/0954-3899/27/3/336
http://iopscience.iop.org/article/10.1088/0034-4885/68/8/R01
http://iopscience.iop.org/0954-3899/46/3/035103
http://iopscience.iop.org/0954-3899/46/3/035103
http://iopscience.iop.org/0954-3899/46/3/035103


Commun. Theor. Phys. 68 (2017) 382–386 Vol. 68, No. 3, September 1, 2017

A Derivation of the Entropy-Based Relativistic Smoothed Particle Hydrodynamics by

Variational Principle

Philipe Mota,1 Wei-Xian Chen (陈渭贤),3 and Wei-Liang Qian (钱卫良)2,3,∗

1Centro Brasileiro de Pesquisas Fisicas, RJ, Brazil

2Escola de Engenharia de Lorena, Universidade de São Paulo, SP, Brazil

3Instituto de F́ısica e Qúımica, Universidade Estadual Paulista Júlio de Mesquita Filho, SP, Brazil

(Received April 17, 2017)

Abstract In this work, a second order smoothed particle hydrodynamics is derived for the study of relativistic heavy
ion collisions. The hydrodynamical equation of motion is formulated in terms of the variational principle. In order to
describe the fluid of high energy density but of low baryon density, the entropy is taken as the base quantity for the
interpolation. The smoothed particle hydrodynamics algorithm employed in this study is of the second order, which
guarantees better particle consistency. Furthermore, it is shown that the variational principle preserves the translational
invariance of the system, and therefore improves the accuracy of the method. A brief discussion on the potential
implications of the model in heavy ion physics as well as in general relativity are also presented.

PACS numbers: 25.75.Ld DOI: 10.1088/0253-6102/68/3/382
Key words: hydrodynamics, SPH method, variational principle

1 Introduction

Hydrodynamics is one of the most venerable theoret-

ical tools which has been playing an important role in

our understanding of nature. Its applications are widely

spread as well as deeply rooted in many distinct areas of

physics. For instance, the hydrodynamic description of

heavy-ion nuclear collisions plays an essential part in the

study of the properties of the hot and dense matter cre-

ated at RHIC and LHC,[1−3] and it is further reinforced

by the onging investigations of fluid/gravity duality.[4−7]

Although the validity and the origin of the hydrodynamic

model have been long under extensive discussions,[8−10]

simulation results[11−15] on azimuthal correlations for var-

ious systems have firmly demonstrated the success of

the approach. The smoothed particle hydrodynamics

(SPH)[16−17] is one of the oldest meshfree methods for the

partial different equation which describes the dynamics of

continuum media. Distinct from any grid-based method

such as the finite element method or the finite difference

method, the SPH makes use of a set of arbitrarily dis-

tributed fluid elements, referred to as particles, to repre-

sent the system. Each particle has a smoothing length, h,

over which their properties are smoothed by a kernel func-

tion. In terms of the kernel function, the contribution of

each particle is weighted according to their distance from

the position in question. Therefore, a physical quantity at

a given spatial point is obtained by summing the relevant

contribution from all the particles lying within the range

of the kernel. The SPH was firstly introduced to study as-

trophysical problems.[16−18] Nowadays, it is widely used to

model fluid motion, as well as solid mechanics.[19]

Despite its wide applications, the original SPH suffers

some inherent problems which lead to low numerical ac-

curacy under certain circumstances. Among others, par-

ticle consistency is one of the notable issues which reflects

the discrepancy between the spatially discretized parti-

cles and the corresponding continuous form of the ker-

nel function. Particle inconsistency demonstrates itself as

the discretized SPH particles to be incapable of properly

reproducing a constant function. It usually results from

the particle approximation process, which is closely associ-

ated with the boundary particles, non-uniformed particle

distribution as well as the smoothing length. The finite

particle method (FPM)[20−21] was proposed by Liu et al .

to improve the particle consistency. The key idea of the

approach is to perform the Taylor series expansion of the

function to be approximated before multiplying both sides

of the equality by the kernel function and integrating over

relevant volume. It was shown that the particle consis-

tency is related to the order of the above Taylor series,

and it is guaranteed independent of the specific form of

the kernel function, neither to the particle distribution.
In implementing the SPH to the partial different equa-

tion, some rules are proposed to symmetrize or asymmet-

ric the terms involving the gradient operator.[22] In the

case of pressure gradient, the term is symmetrized in or-

der to respect Newton’s third law: the pair of forces acting

on the two particles are equal in size but opposite in direc-

tion. Alternatively, it is shown that the above result can

∗E-mail: weiliang qian@gmail.com

c⃝ 2017 Chinese Physical Society and IOP Publishing Ltd

http://www.iopscience.iop.org/ctp http://ctp.itp.ac.cn



No. 3 Communications in Theoretical Physics 383

be obtained naturally, if one derives the hydrodynamic
equation by using the variational principle,[23−25] which
is a consequence that the system conserves linear and
angular momentum. For event by event fluctuating ini-
tial conditions, even though SPH particles are distributed
uniformly at the initial instant, the distribution is likely
to be disturbed as the system evolves in time. There-
fore, FPM formalism is particularly suitable to handle
such physical system. Since the momentum conservation
is important for small systems created in the relativistic
heavy ion collisions, one needs to develop a model which
explicitly preserves the conservation law. Owing to the
complicated form of the FPM, it is not straightforward to
guarantee the momentum conservation by symmetrizing
certain physical quantities. In order to apply the FPM to
relativistic heavy ion collisions, one shall employ the vari-
ational principle to obtain the corresponding equation of
motion. In addition, the system created in the collision is
of significantly high energy density with mostly vanishing
baryon density, therefore the entropy should be chosen as
the base of SPH algorithm. This is the main goal of the
present study. In the following section, we briefly review
the main feature of FPM and discuss its advantage. The
entropy based hydrodynamical equation is derived Sec. 3
by the variational principle. Discussions and conclusions
are given in the last section.

2 The Finite Particle Method
For a physical quantity f(x), the Taylor series expan-

sion gives

f(x) =
∑
n=0

(x− xa)
n
i

n!
(∂n

i f)xa . (1)

Now we multiply both sides by a kernel functionW (x−xa)
and integrate over x to obtain∫

dxf(x)W (x−xa) =
∑
n=0

(∂n
i f)xa

n!

∫
dx(x−xa)

n
i W (x−xa).(2)

If one retains the first term on the right hand side of the

equality∫
dxf(x)W (x− xa) = fa

∫
dxW (x− xa) , (3)

and assumes that
∫
dxW (x − xa) = 1, one restores the

original SPH formula, namely,

fa =

∫
dxf(x)W (x− xa) =

∑
b

νbfb
ρb

W (xb − xa) . (4)

In the last step, one makes use of the particle approxima-

tion.

However, if one applies the following particle approxi-

mation directly to Eq. (3),∫
dxf(x)W (x− xa) →

∑
b

νbfb
ρb

W (xb − xa) , (5)∫
dxW (x− xa) →

∑
b

νb
ρb

W (xb − xa) . (6)

One obtains instead

fa =

∑
b(νbfb/ρb)W (xb − xa)∑
b(νb/ρb)W (xb − xa)

. (7)

It is noted that denominator on the right hand side of

the equation is not exactly “1” in practice, and Eq. (7)

is known as corrective smoothed particle method (CSPM)

in literature, which preserves the zeroth order kernel and

particle consistency.

It is intuitive to generalize the above procedure to

higher order. By retaining the right hand side of Eq. (1)

to the second order, one obtains∫
dxf(x)W (x− xa) = fa

∫
dxW (x− xa)

+ ∂jfa

∫
dx(x− xa)jW (x− xa) . (8)

If one replacesW (x−xa) by [(x− xa)i/|x− xa|]W ′(x−xa)

in the above equation, one has∫
dxf(x)

(x− xa)i
|x− xa|

W ′(x− xa) = fa

∫
dx

(x− xa)i
|x− xa|

W ′(x− xa) + ∂jfa

∫
dx(x− xa)j

(x− xa)i
|x− xa|

W ′(x− xa) . (9)

By implementing particle approximation, Eqs. (8)–(9) cor-
respond to a matrix equation of (D + 1) dimension, with
D being the spatial dimension of the system, as follows[

⟨f⟩a
⟨f⟩a,j

]
=

[
⟨1⟩a ⟨∆xk⟩a
⟨1⟩a,j ⟨∆xk⟩a,j

] [
fa

fa,k

]
, (10)

where

⟨f⟩a ≡
∑
b

νbfb
ρb

Wab , (11)

⟨f⟩a,j ≡
∑
b

νbfb
ρb

(xab)j
|xab|

W ′
ab . (12)

For Eq. (9), in general one may freely replace
W (x − xa) by any basis function, and in particu-
lar, by W ′(x − xa) as done in Ref. [20]. Our
choice of ((x− xa)i/|x− xa|)W ′(x − xa) garantees that

[(x− xa)i/|x− xa|]W ′(x− xa) is an even function as the
kernel W . The above equation can be used to express fa
and ∂jfa in terms of the properties of SPH particles,[

fa

fa,k

]
=

[
⟨1⟩a ⟨∆xk⟩a
⟨1⟩a,j ⟨∆xk⟩a,j

]−1 [ ⟨f⟩a
⟨f⟩a,j

]
. (13)

In one-dimensional case, it gives

fa =
⟨∆x⟩a,x⟨f⟩a − ⟨∆x⟩a⟨f⟩a,x
⟨1⟩a⟨∆x⟩a,x − ⟨1⟩a,x⟨∆x⟩a

, (14)

fa,x =
⟨1⟩a⟨f⟩a,x − ⟨1⟩a,x⟨f⟩a

⟨1⟩a⟨∆x⟩a,x − ⟨1⟩a,x⟨∆x⟩a
. (15)

It is not difficult to see that for a constant function, the
first line of Eq. (15) naturally leads to the constant, while
the second line guarantees a vanishing first order deriva-
tive.



384 Communications in Theoretical Physics Vol. 68

Fig. 1 (Color online) Plots for SPH fit to the superposition of two random Gaussian function by using standard SPH
as well as FPM. The original function is shown in blue dots, while SPH result is denoted by empty red triangles. The
left column is the results for standard SPH, and the right column is those for FPM. The upper panel corresponds to
uniform particle distribution, and lower panel corresponds to random particle distribution. In the calculations, we make
use of 1000 SPH particles with h = 0.002 in all four cases.

In Fig. 1 we show the SPH fit to the superposition of
two random Gaussian functions by using standard SPH
as well as FPM. The upper panel corresponds to uniform
particle distribution, and lower panel corresponds to non-
uniform particle distribution. One sees that even for uni-
form particle distribution, the standard SPH interpolation
cannot properly reproduce the points on the boundary, as
discussed in Ref. [20]. For non-uniform distribution, FPM
is obviously superior to the standard SPH, which is the
case that one frequently encounters in the heavy-ion nu-
clear collisions. The above study is about fit to a given
function, in what follows, we will explore the properties of
the temporal evolution of the system.

3 Hydrodynamic Equation and Temporal
Evolution
The relativistic hydrodynamic equation for ideal fluid

can be obtained by the conservation of energy-momentum
flux,[2]

d

dτ

( (ϵ+ P )

s
γgijv

j
)
− 1

sγ
∂iP = 0 , (16)

where ϵ, P, s are the energy density, pressure and entropy
density in the co-moving frame, vi, γ are the three velocity
and gamma factor of the fluid element, gij is the metric
in Minkowski space.

It is noted that the conservation of the entropy flow is
valid once there is no viscosity, and it is consistent with the
standard form of entropy based SPH formula,[2] namely,

s∗i =
∑
j

νjW (ri − rj ;h) . (17)

In the case of standard SPH, one substitutes the following
symmetrized form for the pressure gradient,[22]

(∂P )i =
∑
j

νjs
∗
i

( Pi

s∗2i
+

Pj

s∗2j

)
∇iW (ri − rj ;h) , (18)

one leads to the following hydrodynamic equation in terms
of degree of freedom of SPH particles.

d

dt

(
νi
Pi + εi

si
γi vi

)
=

∑
j

−νiνj

[ Pi

s∗i
2 +

Pj

s∗j
2

]
∇iW (ri − rj ;h) , (19)

where quantities with a superscript “*” are evaluated in
the laboratory frame and thus evaluated by using the SPH
interpolation or by the equation of state (EoS), they are
related to the corresponding quantity in the co-moving
frame by a gamma factor (eg. s∗i = γsi) due to Lorentz
contraction.

Similarly, in the case of FPM, the pressure gradient on
the right hand side of Eq. (16) can be written as

(∂P )i = Pi,x =
⟨1⟩i⟨P ⟩i,x − ⟨1⟩i,x⟨P ⟩i

⟨1⟩i⟨∆x⟩i,x − ⟨1⟩i,x⟨∆x⟩i
, (20)

where

⟨P ⟩i =
∑
j

νjPj

ρj
Wij , (21)

⟨P ⟩i,x =
∑
j

νjPj

ρj

(xij)

|xij |
W ′

ij . (22)

However, the above hydrodynamic equation does not



No. 3 Communications in Theoretical Physics 385

take into consideration the momentum conservation. Un-
like the case of standard SPH, it is not obvious how to
straightforwardly write down a symmetrized form as in
Eq. (18) to guarantee that the resultant equation of mo-
tion respects the conservation law. By looking closely at
the right hand side of Eq. (19), one observes that it can
be written as ∑

j

f ij , (23)

with

f ij = −νiνj
Pi

s∗i
2 ∇iW (ri − rj ;h) . (24)

Since the kerner function W is an even function, one finds

f ij = −f ji . (25)

In other words, the force excerted on i-th particle by j-
th particle satisfies Newton’s third law. By using the
variational principle, the translational variance of the La-
grangian density implies the momentum conservation or
Newton’s third law. In what follows, we derive the hy-
drodynamic equation by using the variational approach.
Following Ref. [25], the action of the system can be written
as

LSPH({ri, ṙi}) = −
∑
i

(E
γ

)
i
= −

∑
i

νi(ε/s
∗)i , (26)

where Ei is the “rest energy” of the i-th particle.[26] When
applying the variational principle δSSPH = δ

∫
dtLSPH =

0, we note that one has δEi = −PiδVi in the co-moving
frame, Vi = νi/si and δγ = v · δvγ3, which lead to

0 = δSSPH = −
∫

dt
{∑

i

δri ·
d

dt

[
νi

(Pi + εi
si

)
γivi

]
+
∑
i

νiPi

(si∗)2
δsi

∗
}
. (27)

If Eq. (17) were used, one would find∑
i

νiPi

(s∗i )
2
δs∗i =

∑
i,j

( νiPi

(s∗i )
2
νj +

νjPj

(s∗j )
2
νi

)
∇iW (ri − rj ;h)δri , (28)

and consequently Eq. (19).
Now, to calculate the hydrodynamic equation for the

FPM case in a one-dimensional system, we make use of
Eq. (15), namely,

s∗i =
⟨∆x⟩i,x⟨s⟩i − ⟨∆x⟩i⟨s⟩i,x
⟨1⟩i⟨∆x⟩i,x − ⟨1⟩i,x⟨∆x⟩i

, (29)

on the right hand side of Eq. (27). Before carrying out any
explicit calculation, we note that in this case Newton’s
third law is guaranteed since Eq. (29) is translational in-
variant: it remains unchanged if all SPH particles shift the
same amount xi → xi +X. To be specific, for any quan-
tity ai = gi

∑
j tjW

(e)(x i −x j ;h) where W
(e)(x i−x j ;h)

is any even kernel function, it is straightforward to find

δ
(∑

i

ai

)
=

∑
ij

(gitj + gjti)W
(e)′(x i − x j ;h)δxi

≡
∑
j

(f
(e)
ij + f

(e)
ji )δxi . (30)

Similarly, for any quantity bi = hi

∑
j ujW

(o)(x i − x j ;h)

where W (o)(x i−x j ;h) is any odd kernel function, one has

δ
(∑

i

bi

)
=

∑
ij

(hiuj − hjui)W
(o)′(x i − x j ;h)δxi

≡
∑
j

(f
(o)
ij + f

(o)
ji )δxi . (31)

In either case f
(e,o)
ij = −f

(e,o)
ji is satisfied. By a lengthy

but straightforward calculation, one finds the hydrody-
namic equation as follows

d

dt

(
νi
Pi + εi

si
γi vi

)
=

∑
j

f
(n)
ij , (32)

where

f
(n)
ij = −[(l

(n)
i m

(n)
j + (−1)k

(n)

l
(n)
j m

(n)
i )]

×W (n)′(x i − x j ;h) , (33)

with

k
(1,3,6,8)
i = 1, k

(2,4,5,7)
i = 2 ,

l
(1)
i =

Di⟨s⟩i
Bi

, l
(2)
i =

Di⟨∆x⟩i,x
Bi

,

l
(3)
i = −Di⟨s⟩i,x

Bi
, l

(4)
i = −Di⟨∆x⟩i

Bi
,

l
(5)
i = −CiDi⟨∆x⟩i,x

B2
i

, l
(6)
i = −CiDi⟨1⟩i

B2
i

,

l
(7)
i =

CiDi⟨∆x⟩i
B2

i

, l
(8)
i =

CiDi⟨1⟩i,x
B2

i

,

m
(1,3,5,6,7,8)
i =

νi
ρi
, m

(2,4)
i = νi ,

W (1,6)(x i − x j ;h) =
x2
ij

|xij |
W ′(x i − x j ;h) ,

W (2,5)(x i − x j ;h) = W (x i − x j ;h) ,

W (3,8)(x i − x j ;h) = xijW (x i − x j ;h) ,

W (4,7)(x i − x j ;h) =
xij

|xij |
W ′(x i − x j ;h) ,

Bi = ⟨1⟩i⟨∆x⟩i,x − ⟨1⟩i,x⟨∆x⟩i ,
Ci = ⟨∆x⟩i,x⟨s⟩i − ⟨∆x⟩i⟨s⟩i,x ,

Di =
νiPi

(s∗i )
2
. (34)

It is instructive to see how the resulting hydrodynamic

equation reduces to the standard SPH formulae in its

limit. By comparing Eq. (33) with Eq. (24), it is not diffi-

cult to find that Eq. (24) corresponds to the specific term

f
(2)
ij in one-dimensional case when one assumes ⟨1⟩i → 1

and ⟨1⟩i,x → 0. All other terms disappear when one makes

use of the above limit as well as the symmetry of the ker-



386 Communications in Theoretical Physics Vol. 68

nel function, so that either ⟨∆x⟩i → 0 or ⟨∆x⟩i,x → 1
takes place.

4 Concluding Remarks
To summarize, in this work, we make a preliminary

attempt to study the implementation of FPM into the
entropy-based SPH hydrodynamic model. It is argued
that the equation of motion obtained by using variational
principle, though more complicated in its form, is more
suitable to small and/or fluctuating systems where the
conservation law plays a more stringent role in the dy-
namics. We discussed possible implementations and in
particular, derived the hydrodynamic equation of motion
by using variational principle, where the momentum con-
servation of the system is assured. It is shown how the
obtained equation of motion reduces to the standard SPH
form as a limit, which might be instructive to study the
contributions of individual terms when one needs to intro-
duce approximation for practical reasons.

Owing to the observation of the “ridge” effect in two-
particle correlation in relativistic heavy ion collision, the
fluctuating initial conditions play an increasingly impor-
tant role in the hydrodynamical description of nuclear col-
lisions. The AdS/CFT correspondence states the duality

that two apparently distinct physical theories are closely

connected. It provides another insightful viewpoint of hy-

drodynamics as a gradient expansion in the long wave-

length limit. The applications of the SPH algorithm, from

both aspects, strengthen the ongoing studies on heavy-ion

physics as well as on gravitation theory. The improvement

of particle consistency brought by the FPM, therefore, can

be significant in the context of precision and efficiency of

the numerical approach. It is interesting to implement the

obtained equation of motion for realistic collision simula-

tions, which will be carried out in our subsequent study.

Acknowledgments

We are thankful for valuable discussions with Takeshi

Kodama and Yogiro Hama. We gratefully acknowl-

edge the financial support from Fundação de Amparo à

Pesquisa do Estado de São Paulo (FAPESP), Fundação

de Amparo à Pesquisa do Estado de Minas Gerais

(FAPEMIG), Fundação de Amparo à Pesquisa do Estado

do Rio de Janeiro (FAPERJ), Conselho Nacional de De-

senvolvimento Cient́ıfico e Tecnológico (CNPq), and Coor-

denação de Aperfeiçoamento de Pessoal de Nı́vel Superior

(CAPES).

References
[1] P. Romatschke, Int. J. Mod. Phys. E 19 (2010) 1,

arXiv:0902.3663.

[2] Y. Hama, T. Kodama, and O. Socolowski Jr., Braz. J.
Phys. 35 (2005) 24, arXiv:hep-ph/0407264.

[3] R. Derradi de Souza, T. Koide, and T. Kodama, Prog.
Part. Nucl. Phys. 86 (2016) 35, arXiv:1506.03863.

[4] R. Baier, P. Romatschke, D. T. Son, et al., J. High Energy
Phys. 04 (2008) 100, arXiv:0712.2451.

[5] R. Loganayagam, J. High Energy Phys. 05 (2008) 087,
arXiv:0801.3701.

[6] O. Aharony, S. Minwalla, and T. Wiseman, Class. Quant.
Grav. 23 (2006) 2171, arXiv:hep-th/0507219.

[7] V. Lysov and A. Strominger, arXiv:1104.5502 (2011).

[8] P. B. Arnold, J. Lenaghan, G. D. Moore, and L. G.
Yaffe, Phys. Rev. Lett. 94 (2005) 072302, arXiv:nucl-
th/0409068.

[9] M. P. Heller, R. A. Janik, and P. Witaszczyk, Phys. Rev.
Lett. 110 (2013) 211602, arXiv:1302.0697.

[10] P. Romatschke, Eur. Phys. J. C 77 (2017) 21,
arXiv:1609.02820.

[11] J. Takahashi, et al., Phys. Rev. Lett. 103 (2009) 242301,
arXiv:0902.4870.

[12] D. Teaney and L. Yan, Phys. Rev. C 86 (2012) 044908,
arXiv:1206.1905.

[13] R. P. G. Andrade, F. Grassi, Y. Hama, and W. L. Qian,
Phys. Lett. B 712 (2012) 226, arXiv:1008.4612.

[14] P. Bozek and W. Broniowski, Phys. Lett. B 718 (2013)
1557, arXiv:1211.0845.

[15] W. L. Qian, R. Andrade, F. Gardim, F. Grassi, and Y.
Hama, Phys. Rev. C 87 (2013) 014904, arXiv:1207.6415.

[16] L. Lucy, Astrophys. J. 82 (1977) 1013.

[17] R. Gingold and J. Monaghan, Mon. Not. R. Astro. Soc.
181 (1977) 375.

[18] S. Rosswog, New Astron. Rev. 53 (2009) 78,
arXiv:0903.5075.

[19] L. D. Libersky and A. G. Petschek, 395 (1991) 248.

[20] M. B. Liu and G. R. Liu, Appl. Numer. Math. 56 (2006)
19.

[21] M. B. Liu and G. R. Liu, Arch. Computat. Method E. 17
(2010) 25.

[22] J. Monaghan, Annu. Rev. Astron. Astroph. 30 (1992)
543.

[23] R. A. Gingold and J. J. Monaghan, J. Comput. Phys. 46
(1982) 429.

[24] H. T. Elze, Y. Hama, T. Kodama, M. Makler, and
J. Rafelski, J. Phys. G 25 (1999) 1935, arXiv:hep-
ph/9910208.

[25] C. Aguiar, T. Kodama, T. Osada, and Y. Hama, J.Phys.
G 27 (2001) 75, arXiv:hep-ph/0006239.

[26] H. Goldstein, Classical Mechanics, 3 ed. Addison-Wesley,

(2001).