Newton's equation of motion with quadratic drag force and Toda's potential as a solvable one
Daisuke A. Takahashi
nlin.SI physics.class-ph


]Newton's equation of motion with quadratic drag force and Toda's potential as a solvable one





Research and Education Center for Natural Sciences, Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa 223-8521, Japan

daisuke.takahashi@keio.jp

The family of exactly solvable potentials for Newton's equation of motion in the one-dimensional system with quadratic drag force has been determined completely. The determination is based on the implicit inverse-function solution valid for any potential shape, and hence exhaustive. This solvable family includes the exponential potential appearing in the Toda lattice as a special limit. The global solution is constructed by matching the solutions applicable for positive and negative velocity, yielding the piecewise analytic function with a cusp in the third-order derivative, i.e., the jerk. These procedures and features can be regarded as a generalization of Gorder's construction [Phys. Scr. 2015, 90, 085208] to the energy-dissipating damped oscillators. We also derive the asymptotic formulae by solving the matching equation, and prove that the damping of the oscillation amplitude is proportional to t^-1.








Keywords: quadratic drag force, Toda's exponential potential, elliptic functions, piecewise analytic functions 













§ INTRODUCTION

The theory of nonlinear oscillators, while the subject itself is classical and possesses a very long history <cit.>, has been still receiving renewed interest even in recent years, by finding fruitful relationships and applications in the modern scientific topics. For example, the dynamical systems with many degrees of freedom often falls into a single oscillator equation near the bifurcation point, called the normal form <cit.>. The normal form predicts universal scaling laws in various physical quantities which are far from intuitive from the original equation; the examples of scaling behaviors in the nonlinear Schrödinger equation with a spatial inhomogeneity are found in Refs. <cit.> and applied to the superfluidity phenomena in Bose-Einstein condensates. As another important aspect, the nonlinear oscillators play a role as a testbed to verify modern approximation techniques such as the renormalization group perturbation theory <cit.>, homotopy analysis/perturbation methods <cit.>, and so on. 


In the light of the above-mentioned applications, the solvable models with closed-form solutions could be a powerful tool, since the accuracy of the numerical solution can be tested by directly comparing with the exact one. 

In particular, if the system is described by the equations and/or solutions without fundamental mathematical analyticity, the existence of the exact solution becomes further important, because it is not evident whether a number of general mathematical theorems relying on good differentiability are applicable or not. Gorder's work <cit.>, in which the non-differentiable particle motion under the logarithmic potential has been derived, is an illustrative example. As explained below, the problem addressed in this paper also belongs to the same category. 

Here, we clarify the problem which we solve in this paper. Let us consider the following one-dimensional and one-particle equation of motion with the resistive force proportional to the velocity squared (the quadratic drag force):
	
    M ẍ+γ (sgnẋ)ẋ^2+U'(x)=0,

	where M>0 is the mass,  ẍ is the acceleration, ẋ is the velocity, γ>0 is the drag coefficient, and  U(x) is the potential. Such a drag force becomes important when an object in the air has a high speed. Then, the aim of the present paper is phrased as follows: Identify the family of solvable potentials U(x) for Eq. (<ref>) exhaustively, and constructing the explicit solution, elucidate its analytic and asymptotic natures. 

	The sign function sgnẋ in Eq. (<ref>) is in order for the drag force to have an opposite direction to the motion of the particle, and due to this term, the differential equation (<ref>) is not analytic at ẋ=0. Hence, if we are interested in an exact solution, we must solve it for the cases ẋ≷ 0 separately, and join the solutions smoothly. 
	The solution for the uniform gravity U(x)=Mgx can be found in numerous textbooks <cit.>, but we will provide a richer family which are solvable within the elementary and elliptic functions. This family includes Toda's exponential potential appearing in the Toda lattice <cit.> as a special limit. We especially focus on this special case and construct the solutions with their asymptotics. 

The recent relevant studies in similar or modified settings include the pseudo-oscillator equation with the logarithmic potential modeling the path of electrons in a plasma tube <cit.>, the equation without sgnẋ describing the granular materials <cit.> and those with higher-order potentials <cit.>, the projectiles in two dimension <cit.>, the relativistic and Duffing oscillators governed by homotopy methods <cit.>, and the Duffing oscillator with linear damping force <cit.>. In particular, Ref. <cit.> shares a common feature with our present work, since the solution is constructed by matching procedure and the solution becomes singular at the matching point. In our work, we demonstrate that the similar matching method can be applicable even for the energy non-conserving damped system, where the value of the integration constant gradually changes at every matching point, and the cusp appears in the third-order derivative, that is, the jerk. 

 	The organization of the paper is as follows. In Sec. <ref>, we derive a general implicit solution applicable for any potential. In Sec. <ref>, we find a family of solvable potentials including Toda's one.  In Secs. <ref> and <ref>, we construct the solution for Toda's potential. We solve the differential equation for the cases ẋ≷0, and matching them smoothly, we construct the global solution. In Sec. <ref>, we provide an asymptotic behavior of the solution. Section <ref> is devoted to the discussion, summary, and future outlook.



§ IMPLICIT SOLUTION FOR GENERAL U(X)

We first derive an implicit solution to Eq. (<ref>) for general U(x). 
Since Eq. (<ref>) does not include t explicitly, if x(t) is a solution,  x(t-t_0) is also a solution. Therefore, if we rewrite the equation with respect to not x(t) but its inverse function t(x),  it contains only  t'(x) and t”(x). Defining s(x):=t'(x), and using the formulas for inverse functions ẋ=1/s and ẍ=-s'/s^3, the resultant equation is
	
    -Ms'+γ |s|+s^3U'(x)=0.

	Thus, the second-order equation is reduced to the first-order one. This equation can be linearized by substituting s=± f^-1/2 for s ≷ 0, and we obtain the solution in an implicit form
	
    t= √(M/2)∫±e^±κ x/2dx/√(C-∫e^±κ xU'(x)dx),  κ:=2γ/M.    (ẋ≷0)

	If κ=0, the physical interpretation of the constant C becomes the energy E, and the solution reduces to the well-known form t=√(M/2)∫[E-U(x)]^-1/2dx (e.g., chapter III of Ref. <cit.>). 

	The solution (<ref>), though it is applicable for any U(x), is not actually so convenient, unless the solution can be obtained in an explicit form “x(t)=…”. The reason is as follows.  In the energy-conserving system (i.e., κ=0), even if the integration cannot be performed explicitly for general U(x), it is still useful to identify the region where the particle motion is possible for a given E. Furthermore, if the solution is given by a periodic oscillation, the global solution for fixed E can be constructed by repeating the copy-and-paste of the one-period solution. On the other hand, in the present system, we need to construct a global solution by connecting the solutions for ẋ≷0, where the value of the constant C changes at every junction point ẋ=0 because of the damping by drag force. Hence, the solution cannot be obtained by a naive repetition of the one-period solution.  Therefore, finding a potential such that the explicit closed-form solution is available for arbitrary C is essential.
	


§ FINDING A SOLVABLE POTENTIAL

	Let us now find a solvable potential such that the integration in Eq. (<ref>) can be performed explicitly. 


	Introducing the new variables z=e^±κ x for the cases ẋ≷0, Eq. (<ref>) is rewritten as
	
    t=√(M/2κ)∫^exp(κ x)dz/√(z[ Cκ -∫ u(z)dz ])    (ẋ>0), 
    √(M/2κ)∫^exp(-κ x)dz/√(z[ Cκ+∫ u(z^-1)dz ])    (ẋ<0),

	where u(z):= U'(ln z/κ).  We now seek a potential such that Eq. (<ref>) can be calculated within the elementary and elliptic functions. It imposes the condition that both z[ Cκ -∫ u(z)dz ] and z[ Cκ+∫ u(z^-1)dz ] are polynomials of order ≤ 4, and it fixes the form of the potential as u(z)=a+bz^2+cz^-2. The corresponding U(x) is given by
	
    U(x)=A/2κ( 2κ xtanh2κ d+cosh 2κ(x-d)/cosh 2κ d-1 )

	with A ∈ℝ and d∈ℝ or ℝ+iπ/4κ. Here, we chose U(0)=U'(0)=0 without loss of generality. In this normalization,  u(z) is written as
	
    u(z)=A( z_0^-2z^2-z_0^2z^-2/z_0^2+z_0^-2+z_0^2-z_0^-2/z_0^2+z_0^-2)

	with d=ln z_0/κ. 

	If d ∈ℝ and A>0, the potential (<ref>) is non-negative everywhere and the particle motion is always bounded. If d ∈ℝ+iπ/4κ, Eq. (<ref>) has local extrema at x=0 and 2d, and the motion may not be bounded depending on the initial condition. If we set d → +∞, it reduces to
	
    U(x)=A( x+e^-2κ x-1/2κ),

	which is just the famous potential appearing in the Toda lattice <cit.>. The emergence of Toda's exponential potential in the present context is a little unexpected, since we have made no soliton-theoretical consideration here. In this case the corresponding u(z) is given by u(z)=A(1-z^-2). Henceforth, we concentrate on this special potential and investigate its solutions in detail.



§ SOLUTIONS OF THE DIFFERENTIAL EQUATIONS FOR Ẋ≷0

	In this and the next section, we construct the solution for the potential (<ref>). 
	Among the solvable family in Eq. (<ref>), this potential is the easiest one in the sense that the solution for ẋ>0 becomes an elementary function. 
 
	This section provides the list of the solutions to the differential equations for ẋ≷0, which will become the “pieces” of the global solution given in Sec. <ref>. Henceforth we follow the convention by Abramowitz and Stegun for the notation of the elliptic integrals and functions. 

	For u(z)=A(1-z^-2), Eq. (<ref>) reduces to
	
    t = √(M/2κ A)∫^exp(κ x)dz/√(2c_+ z-1-z^2)    (ẋ>0), 
    √(3M/2 κ A)∫^exp(-κ x)dz/√(z(c_-+3z-z^3))    (ẋ<0),

	where c_+=κ C/2A and c_-=3 κ C /A. There are one real-valued solution for ẋ>0 and three for ẋ<0 depending on the values of constants c_±. Below we summarize it.  

	I. The case ẋ>0. The real-valued solution exists when c_+ ≥ 1. Writing c_+=coshκ x_0 with x_0≥0, the solution is given by
	
    x_I(t,x_0):=1/κln[ coshκ x_0-(sinhκ x_0)cosω_I t ], 
       ω_I=√(2κ A/M).
 
	II. The case ẋ<0 and -2 ≤ c_-≤ 0. Parametrizing c_-=-2cos3φ, φ∈ [0,π/6], the polynomial is factorized as z(c_-+3z-z^3)=-z∏_n=0,1,2[z+2cos( φ-2nπ/3 )]. The resultant solution is then given by
	
    x_II(t;φ):=1/κln[ sn^2(ω_IIt/2|m_II)/2sin(π/6+φ)+cn^2(ω_IIt/2|m_II)/2sin(π/6-φ)], 
       m_II=m_II(φ):=sin2φ/cos(2φ-π/6),
       ω_II=ω_II(φ):=ω_I√(2cos(2φ-π/6)/√(3)).
  
	III. The case ẋ<0 and 0 ≤ c_-≤ 2. Parametrization is the same as II, but φ∈ [π/6,π/3]. The solution is then given by
	
    x_III(t;φ):=1/κln[ ns^2(ω_IIIt/2|m_III)/2sin(π/6+φ)-cs^2(ω_IIIt/2|m_III)/2sin(π/6-φ)],
       m_III=m_III(φ):=cos(2φ-π/6)/sin2φ,
       ω_III=ω_III(φ):=ω_I√(2sin2φ/√(3)).

	IV. The case ẋ<0 and 2 ≤ c_-. Parametrizing c_-=2cosh3φ, φ≥0, the solution is  
	
    x_IV(t;φ):=
       1/κln[ 1/2coshφ( 1+√(3(2cosh2φ+1)/2cosh2φ-1)1-cn(ω_IVt|m_IV)/1+cn(ω_IVt|m_IV)) ],
       m_IV=m_IV(φ):=1/2-√(3)/2cosh2φ/√(1+2cosh4φ),
       ω_IV=ω_IV(φ):=ω_I[1+2cosh4φ/3]^1/4.

	When c_-<-2, no real solution exists. 

	Note that the solutions  x_II,x_III, and x_IV reduce to the same expression if the modulus m=k^2 ∈ [0,1] is allowed to take more general values. For example,  x_III reduces to x_II using ns(z+iK'|m)=sn(√(m)z|m^-1). The complete classification of real-valued elliptic integrals expressed only using real parameters and the modulus ∈ [0,1] is given in Ref. <cit.>. (Setting z=-t, use 253.00 for  x_II and x_III, and 259.00 for  x_IV.)



§ GLOBAL SOLUTION

	We now provide the global solution for the potential (<ref>) applicable for all time t≥0, which can be constructed by smoothly matching the solutions x_I, x_II, x_III, and x_IV at the points ẋ(t)=0. For simplicity, here we consider the specific initial condition x(0)=-x_0<0 and ẋ(0)=0. In this case the solution is given by x_I→ x_II→ x_I→ x_II→…. When x(0)>0 and ẋ(0)=0, x_II comes first and the rest of the motion is the same. Note that x_III and x_IV are necessary only for the initial condition with very large ẋ(0)<0 (or very large x(0)<0 with not-so-large ẋ(0)<0), where the high-speed particle climbs the exponentially increasing potential slope.  They are used only once; after starting from x_III (or x_IV), the remaining motion is described by x_III (or x_IV) → x_I→ x_II→ x_I→ x_II→…. 








	Here we show the result. The solution satisfying the initial condition x(0)=-x_0<0, ẋ(0)=0 is given by
	 
    x(t)=
    			x_I( t-t_n,x_n )     t_n ≤ t ≤ t_n+T_I, 
    
    			x_II( t-t_n-T_I,φ_n )     t_n+T_I≤ t ≤ t_n+1,

	where T_I,  t_n, x_n, and φ_n are defined as follows. Let us write
	
    T_I=π/ω_I,   T_II(φ)=2K(m_II(φ))/ω_II(φ),

	where ω_I, ω_II(φ), and m_II(φ) are introduced in the previous section, and  K(m) is the complete elliptic integral of the first kind. Let us define φ_0 ∈ [0,π/6) by the relation x_0=-1/κln[2sin(π/6-φ_0)], and let us define φ_1, φ_2,… by the recurrence relation
	
    sin( π6-φ_n+1)=1/4sin( π6+φ_n ),

	and write x_n:=-1/κln[2sin(π/6-φ_n)]. Equation (<ref>) arises from the matching condition between x_I and x_II.  Then,  t_n's are given by 
	
    t_0=0,   t_n:=nT_I+∑_j=0^n-1T_II(φ_j),   (n=1,2,3,…).

	Here we briefly sketch the derivation of Eq. (<ref>) by the matching condition. Setting t=t_n+1 in Eq. (<ref>), we obtain x_I(0,x_n+1)=x_II(T_II(φ_n),φ_n)  ↔ x_n+1=-1/κ[2sin(π/6+φ_n)]. On the other hand, if we set t=t_n+1+T_I, we have x_I(T_I,x_n+1)=x_II(0,φ_n+1)  ↔  -x_n+1=-1/κ[2sin(π/6-φ_n+1)]. Eliminating  x_n+1 from these two equations, we obtain (<ref>). 

	The plot of the solution is given in Fig. <ref>.  
	Note that  x(t), ẋ(t), ẍ(t) and ⃛x(t) are continuous, while the higher-order derivatives d^n/dt^nx(t)  (n≥ 4) are generally discontinuous at the points ẋ=0, originating from the non-differentiable sgnẋ factor in Eq. (<ref>). We can indeed observe the cusp of the jerk ⃛x(t) in Fig <ref>.
	



§ ASYMPTOTICS AND ENVELOPES

	
	Here, we determine the asymptotic behavior of the solution (<ref>). The solution of the matching condition, which is given by the recurrence relation (<ref>), is well approximated by
	
    sin( π/6-φ_n ) ≃1/2(1-3[1-2sin( π/6-φ_0 )]/3+4n[1-2sin( π/6-φ_0 )]).

	It follows from that δ_n:=1/2-sin( π/6-φ_n ) satisfies δ_n+1-δ_n ≃8/3δ_n^2. Using Eq. (<ref>), the approximate expression for  x_n=-1/κln[2sin(π/6-φ_n)] is
	
    x_n ≃3x_0/3+4nκ x_0,

	if the initial x_0 is not so large (κ x_0 ≲ 1). Within the same approximation,  T_II(φ_j) is estimated to be
	
    T_II(φ_n)≃ T_I[ 1+5/12( 3κ x_0/3+4nκ x_0)^2 ].

	Using this, the summation of t_n in Eq. (<ref>) can be carried out using the polygamma function. For large n, it reduces to
	
    t_n≃ 2nT_I[ 1+5κ^2 x_0^2/8(3+4nκ x_0)].

	Solving it with respect to n, we get  n≃t_n/2T_I( 1-5κ x_0/16T_I/t_n), and using this, we can eliminate n from Eq. (<ref>). Thus we arrive at the expression for the lower envelope




	
    x_envelope(t)= -3T_Ix_0/3T_I+2κ x_0 t.

	The upper envelope is given by -x_envelope(t-T_I).

	Next, let us derive the asymptotic curve. Within the cosine-curve approximation, the n-th oscillation is given by ∝cos2π(t-t_∞)/t_n-t_n-1. Up to the present accuracy, we simply obtain t_n-t_n-1≃ 2T_I. The phase shift t_∞ is determined by imposing that the argument of the cosine becomes 2nπ at t=t_n, and we have t_∞=5T_Iκ x_0/16. Summarizing, the asymptotic curve is
	
    x_asymptotic(t)=x_envelope(t) cos[ π( t/T_I-5κ x_0/16) ].

	The verification of these envelopes and the asymptotic curve is shown in Fig. <ref>. 

	The damping by t^-1 shown in Eq. (<ref>) is very slow compared to the system with linear drag force, where the decay occurs in an exponential way. This comes from the fact that the quadratic drag force is very weak when the particle moves slowly. 





§ DISCUSSION, SUMMARY, AND FUTURE OUTLOOK

	The main findings of the present paper can be summarized by the following (a) and (b):


  * Determination of the family of solvable potentials (<ref>) for the equation of motion with the quadratic drag term (<ref>).

  * Generalization of the matching method to construct the piecewise analytic solution in Ref. <cit.> to the damped oscillators.

The result (a) can be accomplished using the general implicit solution (<ref>) and (<ref>) which are valid for arbitrary potentials. Moreover, the family includes the Toda potential (<ref>) as a limiting case. The global solution for this potential is constructed by matching the analytic solutions valid for each ẋ≷0, and the resultant solution has a cusp in the third-order derivative. Such a matched solution with piecewise analytic property is the same as Ref. <cit.>. The main new difference of our result is stated in (b). Namely, the value of the integration constant C in Eq. (<ref>) changes at every junction point ẋ=0, where the generalized matching procedure well describes the damped oscillation in the dissipative system. This new situation increases the importance of the use of the solvable potential. Furthermore, we have also provided the asymptotic curves of the solutions by solving the matching condition, and have shown that the damping occurs by the rate proportional to t^-1 [Eqs. (<ref>) and (<ref>)]. 

Finally, we provide a few prospects. (i) In this paper, the exhaustive determination of solvable potentials has been made within the elementary and elliptic functions. If we want to extend the argument using more general class of functions such as hypergeometric or hyperelliptic ones, Eq. (<ref>) will still remain to be a starting point due to its generality. (ii) The investigation of the solutions for other solvable potentials in Eq. (<ref>) is left as a future task. Since they have local extrema and solutions for both signs of velocity becomes elliptic,  we expect more complicated classification based on their stability and boundedness. Verifying the damping behavior of the solutions [Eqs. (<ref>) and (<ref>)] for these more general potentials will be also important. For example, if we set d=0 in Eq. (<ref>), we get a potential with symmetric shape U(x)=A/2κ(cosh 2κ x-1), which will be more suitable from the viewpoint of experimental realization and investigation of damping dynamics. One of the most well-known system described by Eq. (<ref>) is the classical ball with air's drag. So, if we construct a bowl with the cosh-like shape or the Toda-type shape [Eq. (<ref>)], and consider the motion of the ball on this bowl in the atmosphere, the mathematical result in this paper will be verified experimentally.  (iii) Whether the appearance of Toda's potential has some implication to the soliton theory might also be worth considering. In general, soliton equations possess the inhomogeneous and higher-order generalizations. The former arises from the non-isospectral problem and used to describe relaxation phenomena in a nonuniform medium <cit.>, and the latter is made from the family of the Lax pairs in the hierarchy <cit.>, which were recently used as an effective model of density-modulated quantum condensates <cit.>. If we could invent a new dissipative effect for these equations with preserved integrability by extending the idea of the present paper, it would be useful to study damping phenomena of nonlinear wave packets.

This work is supported by the Ministry of Education, Culture, Sports, Science (MEXT)-Supported Program for the Strategic Research Foundation at Private Universities “Topological Science” (Grant No. S1511006). 



iopart-num



10
url<#>1#1urlprefixURL 



LambDynamics
Lamb H 1923 Dynamics 2nd ed (Cambridge: Cambridge University Press)

WhittakerTreatiseOn
Whittaker E T 1937 A Treatise On the Analytical Dynamics of Particles and
  Rigid Bodies 4th ed (Cambridge: Cambridge University Press)

GuckenheimerHolmes
Guckenheimer J and Holmes P 1983 Nonlinear oscillations, dynamical
  systems, and bifurcations of vector fields (Berlin: Springer)

PhysRevE.55.2835
Hakim V 1997 Phys. Rev. E 55(3) 2835–2845

HUEPE2000126
Huepe C and Brachet M E 2000 Physica D 140 126–140

PHAM2002127
Pham C T and Brachet M 2002 Physica D 163 127–149

TAKAHASHI20121589
Takahashi D A 2012 Physica D 241 1589–1602

PhysRevLett.105.035302
Kato Y and Watabe S 2010 Phys. Rev. Lett. 105(3) 035302

PhysRevE.54.376
Chen L Y, Goldenfeld N and Oono Y 1996 Phys. Rev. E 54(1)
  376–394

doi:10.1143/PTP.94.503
Kunihiro T 1995 Prog. Theor. Phys. 94 503–514

doi:10.1142/S0217979206033796
He J H 2006 Int. J. Mod. Phys. B 20 1141–1199

LiaoHAM
Liao S 2013 Advances in the Homotopy Analysis Method (Singapore: World
  Scientific)

1402-4896-90-8-085208
Gorder R A V 2015 Phys. Scr. 90 085208

Todalattice
Toda M 1989 Theory of Nonlinear Lattices 2nd ed (Berlin: Springer)

1402-4896-89-10-105205
Gadella M and Lara L P 2014 Phys. Scr. 89 105205

PhysRevE.51.2538
Linz S J and Hänggi P 1995 Phys. Rev. E 51(3) 2538–2542

0143-0807-16-2-004
Linz S J 1995 Eur. J. Phys. 16 67–72

1402-4896-85-4-045006
Lai S K and Chow K W 2012 Phys. Scr. 85 045006

doi:10.1119/1.10812
Parker G W 1977 Am. J. Phys. 45 606–610

1751-8121-40-29-015
Yabushita K, Yamashita M and Tsuboi K 2007 J. Phys. A: Math. Theor.
  40 8403

1402-4896-77-2-025004
Beléndez A, Pascual C, Fernández E, Neipp C and Beléndez T 2008
  Phys. Scr. 77 025004

0143-0807-36-6-065020
Johannessen K 2015 Eur. J. Phys. 36 065020

LLMechanics
Landau L D and Lifshitz E M 1976 Mechanics 3rd ed (Oxford:
  Butterworth-Heinemann)

ByrdFriedman
Byrd P F and Friedman M D 1971 Handbook of Elliptic Integrals for
  Engineers and Scientists 2nd ed (Berlin: Springer)

PhysRevLett.37.693
Chen H H and Liu C S 1976 Phys. Rev. Lett. 37(11) 693–697

doi:10.1143/JPSJ.41.2141
Hirota R and Satsuma J 1976 J. Phys. Soc. Jpn. 41 2141–2142

GUPTA1979420
Gupta M 1979 Phys. Lett. A 72 420–422

FaddeevTakhtajan
Faddeev L D and Takhtajan L A 1987 Hamiltonian Methods in the Theory of
  Solitons (Berlin: Springer)

CORREA20092522
Correa F, Dunne G V and Plyushchay M S 2009 Ann. Phys. 324
  2522–2547

TAKAHASHI2012632
Takahashi D A, Tsuchiya S, Yoshii R and Nitta M 2012 Phys. Lett. B 718 632–637

1402-4896-90-4-045205
Liu D Y, Tian B, Sun W R and Wang Y P 2015 Phys. Scr. 90 045205

PhysRevE.93.062224
Takahashi D A 2016 Phys. Rev. E 93(6) 062224




