Exact solution of Smoluchowski's equation for reorientational motion in Maier-Saupe potential
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t Exact solution of Smoluchowski’s equation forreorientational motion in Maier-Saupepotential
A.E. Sitnitsky,
Institute of Biochemistry and Biophysics, P.O.B. 30, Kazan 420111, Russia.e-mail: [email protected]
Abstract
The analytic treatment of the non-inertial rotational diffusion equation, i.e., of theSmoluchowski’s one (SE), in a symmetric genuinely double-well Maier-Saupe uni-axial potential of mean torque is considered. Such potential may find applicationsto reorientations of the fragments of structure in polymers and proteins. We obtainthe exact solution of SE via the confluent Heun’s function. The solution is uniformlyvalid for any barrier height. We apply the obtained solution to the calculation ofthe mean first passage time and the longitudinal correlation time and obtain theirprecise dependence on the barrier height. In the intermediate to high barrier (lowtemperature) region the results of our approach are in full agreement with thoseof the approach developed by Coffey, Kalmykov, D´ejardin and their coauthors. Inthe low barrier (high temperature) region our results noticeably distinguish fromthe predictions of the literature formula and give appreciably greater values for thetransition rates from the potential well. The reason is that the above mentionedformula is obtained in the stationary limit. We conclude that for very small bar-rier heights the transient dynamics plays a crucial role and has to be taken intoaccount explicitly. When this requirement is satisfied (as, e.g, at the calculation ofthe longitudinal correlation time) we obtain absolute identity of our results withthe literature formula in the whole range of barrier heights. The drawbacks of ourapproach are its applicability only to the symmetric potential and its inability toyield an analytical expression for the smallest non-vanishing eigenvalue.
Key words: rotational motion, diffusion, confluent Heun’s function.
Email address: [email protected] (A.E. Sitnitsky).
Preprint submitted to Physica A 9 June 2018
Introduction
Rotational reorientations are a particular type of motion whose utmost impor-tance for applications can hardly be underestimated. They are ubiquitous inphysics-chemical studies of liquid crystals, polymers, proteins, lipids and manyother kinds of the stuff. Rotational reorientations are investigated experimen-tally with the help of NMR, dielectric relaxation spectroscopy, fluorescencedepolarization, etc (see, e.g., [1] and refs. therein). The main theoretical toolfor their investigation is the Smoluchowski’s equation (SE), i.e., non-inertialrotational diffusion equation for a rigid body in an external potential of meantorque. It is an approximation to the more general master equation (ME).The latter is one of the basic tools in non-equilibrium statistical mechanicsand physical kinetics [2], [3], [4]. ME is that for the time evolution of theprobability density and consequently it expresses the fundamental principleof kinetic balance. Continuous ME is an integro-differential equation and gen-erally it is rather difficult for analytical treatment. ME is the starting pointfor all models of molecular reorientation [5]. In each particular case physicalintuition and/or first principles calculations have to be used in order to formu-late an explicit expression for the transition probability which determines theentire process. However, in only very few examples the transition probabilityallows an analytic solution of ME. This difficulty is the reason why (insteadof looking for an analytic solution) one usually tries to find approximations tothe original ME. The well-known Kramers-Moyal expansion, e.g., transformsME into a partial differential equation of infinite order. The transition fromME (by Taylor series expansion of the jump probabilities and the probabilitydensity for infinitely small jump steps) results in the Fokker-Planck equationwhich is a partial differential one [2], [3], [6]. For the case when the inertialeffects are negligible (overdamped limit) the latter is reduced to SE. Thereis vast literature on their study and applications [2], [3], [6], [7], [8], [9], [10],[11], [12], [13], [14], [15], [16], [17], [18], [20], [21], [22], [23], [24], [25], [26], [27],[28], [29]. The applications of them to the dielectric spectroscopy [30], [31],[27], fluorescence depolarization [16] and NMR relaxation of liquid crystals(see [32], [33], [34] and refs. therein) have been thoroughly explored. Also,the extension of the theory from ordinary diffusion to the fractional one isintensively studied (see [26], and refs. therein).As mentioned above, SE is only an approximation to ME valid for the caseof infinitely small jump steps in a space variable. According to [5] all modelsfor molecular motion are divided into a jump process and a diffusion process.The jumps can take place between discrete set of accessible places or in thecontinuous region of accessible places. In the latter case the diffusion limitis that of infinitely small jump steps. There is a fundamental result that ”adiffusion process always be approximated by a jump process, not the reverse”[3]. Thus the diffusion model ([8], [10], [11], [16], [21], [22]) can be inferred from2he jump model. In this regard SE is subordinate to ME. One has to resortto SE because in the general case its treatment is easier than that of ME. SEarises in the theory of ferromagnetism (where it is called Brown’s equation)[37] and that of molecular rotational motion in a uniaxial potential [28], [26](occurring for dynamics of liquid crystals, dielectric relaxation, etc.). In itsturn, a particular case of SE in a uniaxial double-well potential of the meantorque can be analyzed with the help of the effective potential comprisingthe Maier-Saupe one as an ingredient. The latter is widely used in the theoryof rotational reorientations of nematic liquid crystals [32], [33], [34], theoryof ferromagnetism [35], [36], [37] and in the theory of dielectric relaxation ofrod-like molecules [26], [28].The most powerful approach to the problem was developed by Coffey, Kalmykov,D´ejardin and their coauthors in the above mentioned papers [7], [17], [26], [27],[28], [29]. For the sake of brevity we further call it as CKD. The approach isbased on the expansion of the probability distribution function as a series ofspherical harmonics. This method is well suited for the potentials of meantorque that can be expanded in terms of spherical harmonics. It results in aninfinite hierarchy of differential-recurrence relations for the the moments (theexpectation values of the spherical harmonics). CKD is an exact approach thatuses no approximations and imposes no physical limitations. For instance, theexact solution in the form of the Green function in the frequency domain wasgiven as a continued fraction [38]. Moreover, the correlation time, mean firstpassage time, etc. may be rendered exactly by writing the continued fractionsin integral form. Also, the finite integral representations of the various relax-ation times may be calculated directly by quadratures [39], [40], [41]. However,CKD still leaves room for complementary development that uses the expan-sion of the solution over eigenfunctions of Smoluchowski’s operator rather thanthat over Legendre polynomials.In the present paper we deal with one-dimentional motion along the polarangle ψ coordinate. Such situation can arise in two ways. First, the system isinvariant relative the rotations along the azimuthal angle ϕ coordinate. Thiscase is familiar in the theory of nematic liquid crystals, that of ferromagnetism,etc. and is well described by the ordinary Maier-Saupe uniaxial potential ofmean torque having the minimums at ψ = 0 and ψ = π . Many well knownresults exist for this case in the literature (see above cited papers). The nor-malization of the probability distribution function for this case is given below(1). Second, the system undergoes reorientations along the polar angle ψ co-ordinate at fixed azimuthal angle ϕ = ϕ = const . Such case can take placefor the motion (in molecular frame) of fragments of structure in polymers orin proteins. For instance, the motion of Ω-loops in proteins (see, e.g., [42], Sec.III A in [43] and refs. therein) may be considered as a reorientation betweenstable conformations one of which is ”open” state and the other is ”closed”one. Such motion has important functional role because the Ω-loop usually3lays the role of a lid that regulates the penetration of a ligand into proteininterior or the access of a substrate into enzyme active site as, e.g., in thecase of triosephosphate isomerase [42]. These reorientations between stableconformations may be considered as a transition between the wells of the cor-responding potential over its barrier. It should be stressed that the minimumsof the potential in this case do not necessarily lie at ψ = 0 and ψ = π andwe need a genuinely double-well potential. A convenient way to obtain suchpotential is to add a logarithmic contribution to the Maier-Saupe one. Suchcontribution is suggested by the potential devised in the paper of Pastor andSzabo [44]. There may be different suitable functions under the logarithm butwe show that the one used in the present paper is the mostly convenient onebecause it makes it possible to obtain the representation of the eigenfunctionsof the Smoluchowski’s operator via the known and rather well studied (up tothe point that it is tabulated in Maple) special function, namely the confluentHeun’s function (CHF). The resulting double-well potential for reorientationalmotion is a model typical one analogous to the famous − ax + bx for trans-lational motion. Unfortunately, there is a tiny problem. The normalization ofthe probability distribution function for this case is14 π π Z dϕ δ ( ϕ − ϕ ) π Z dψ sin ψ f ( ψ, t ) = 14 π π Z dψ sin ψ f ( ψ, t ) = 1In all other respects the description of the motion along the polar angle ψ coordinate remains absolutely the same as in the first case. Thus, the resultsfor two realistic cases can not be directly compared with one another becauseit is necessary to keep in mind some ”transformation” constant. We find itawkward to suggest a model for the second case that can not be directlycompared with the previous literature data for the first case. For this reasonwe construct a combined model potential. We use the double well Maier-Saupepotential of the second case for the polar angle ψ but at the same time usethe isotropy over the azimuthal angle ϕ of the first one (i.e., make use of thenormalization (1)). The resulting uniaxial potential of mean torque may seemunphysical but it has the advantage that it provides direct comparison of ourresults with the known literature ones especially with those of CKD approach.We show that the approach developed in the present paper is valid for suchhypothetical potential and hence it remains robust for a realistic situation ofthe second case.Thus, it seems interesting to obtain the analytic solution of the problem mak-ing use the expansion over eigenfunctions of Smoluchowski’s operator for thedouble well Maier-Saupe potential that makes it possible and to compare theresults of CKD with it. In the present paper we show that for the case ofsymmetric potential the problem under consideration appears to be amenableto stringent analytic treatment. We obtain the exact solution for the case of4E for reorientational motion in a symmetric double-well Maier-Saupe uniax-ial potential of mean torque via CHF. Our solution is uniformly valid for anybarrier height. The CHF is a known and by now well described special func-tion which is a solution of the confluent Heun’s equation [45], [46], [47], [48].We apply the obtained solution to the calculation of the mean first passagetime (MFPT) and the longitudinal correlation time and obtain their precisedependence on the barrier height.The paper is organized as follows. In Sec. 2 the problem under study is for-mulated. In Sec. 3 the solution of SE is presented. In Sec. 4 the probabilitydistribution function is obtained. In Sec. 5 the general result is exemplified bythe calculation of the escape rate from a well in the double-well Maier-Saupeuniaxial potential of mean torque. Also the longitudinal correlation time iscalculated. In Sec. 6 the results are discussed and the conclusions are summa-rized. In Appendix the numerical calculations with CHF are considered. In the spherical frame Ξ ≡ { ψ, ϕ } (where 0 ≤ ψ ≤ π is the polar angle and0 ≤ ϕ ≤ π is the azimuthal one) we introduce the conditional probability P (Ξ , t ; Ξ ,
0) of finding the probe at orientation Ξ at time t , if the orientationwas Ξ at time zero and the equilibrium distribution function P eq (Ξ). Thelatter represents the equilibrium probability of finding the probe at orienta-tion Ξ and is connected to the anisotropic potential of mean torque U (Ξ)through the Boltzmann distribution. In many practical situations the distri-bution function P (Ξ , t ; Ξ ,
0) and P eq (Ξ) characterizing a system of interestusually depend only on the polar angle ψ (that between the axis of the probeand the z axis of the chosen frame). The most notable example is nematicliquid crystals in a uniaxial phase for which a rotation about the director,assumed to be the z laboratory axis, should leave the system invariant. Suchsituations arise for axially symmetric potentials when there are no dynamicalcoupling between the longitudinal and the transverse modes of motion. In thiscase the longitudinal modes are governed by the single state variable (polarangle ψ that is called colatitude [37]) while the azimuthal angle gives rise onlyto a steady precession round the axis z . As is stressed in [37] in the theory offerromagnetism the exact Fokker-Planck equation in the single variable (i.e.,SE) follows directly from the axial symmetry of the potential. In the theoryof non-inertial rotational diffusion for a rigid body the requirement of strongdamping is necessary besides the above one [37].For the sake of brevity we further denote such distribution functions as P (Ξ , t ; Ξ ,
0) = P ( ψ, t ; ψ , ≡ f ( ψ, t )5nd P eq (Ξ) = P eq ( ψ ) ≡ f eq ( ψ ) = f ( ψ, t → ∞ ). The latter is defined by thepotential of mean torque V ( ψ ) originating from the long range order of thesystem f eq ( ψ ) = const exp [ − V ( ψ ) / ( k B T )]The distribution function f ( ψ, t ) must be a solution of ME or approximatelyof the corresponding SE. The latter is the equation for the time evolution ofthe distribution function f ( ψ, t ). This function must be normalized so that itsintegral over the whole space gives14 π π Z dϕ π Z dψ sin ψ f ( ψ, t ) = 12 π Z dψ sin ψ f ( ψ, t ) = 1 (1)The normalized initial condition takes the form f ( ψ,
0) = 2sin ψ δ ( ψ − ψ ) (2)For a uniaxial potential of mean torque V ( ψ ) SE under consideration is [28],[26] 2 τ ∂f ( ψ, t ) ∂t = β sin ψ ∂∂ψ " sin ψf ( ψ, t ) ∂V ( ψ ) ∂ψ +1sin ψ ∂∂ψ " sin ψ ∂f ( ψ, t ) ∂ψ (3)where τ is the characteristic relaxation time for isotropic non-inertial rota-tional diffusion (e.g., τ D the Debye one for the theory of dielectric relaxationor τ N the N´eel one for the theory of ferromagnetism) and β = 1 / ( k B T ). Weintroduce a new variable x = cos ψ (4)Then the equation takes the form ∂f ( x, t ) ∂t = 12 τ ∂∂x " (1 − x ) ∂f ( x, t ) ∂x + βf ( x, t ) V ′ ( x ) ! (5)where the dash means the derivative over variable x .6n the stationary limit t → ∞ the probability distribution function f ( x, t )must tend to its equilibrium value f eq ( x ). Thus we have ∂f ( x, t ) ∂t = 0and f ( x, t → ∞ ) → f eq ( x ) = const exp − V ( x ) k B T ! (6)Further we consider the effective double-well potential with logarithmic con-tribution − β ln(1 − x ). This form is suggested by the potential devised in[44] V ( x ) = U ( x ) − β ln(1 − x ) (7)where the ingredient U ( x ) is responsible for the barrier. The Maier-Saupecontribution is most widely used in the literature U ( x ) = U MS ( x ) = b (1 − x ) (8)where the parameter b defines the barrier height. We introduce the dimension-less parameter α = βb U MS ( x ). Equation (5) can be solved by separation of variables f ( x, t ) = η ( x ) χ ( t ) (10)Denoting the separation constant as − ρ ( ρ >
0) we obtain χ ( t ) = exp ( − ρt ) (11)7 .0 0.5 1.0 1.5 2.0 2.5 3.0 Ψ H radians L H Ψ L(cid:144)H k B T L b (cid:144)H k B T L = (cid:144)H k B T L = (cid:144)H k B T L = (cid:144)H k B T L = (cid:144)H k B T L = (cid:144)H k B T L = Fig. 1. The symmetric double-well Maier-Saupe uniaxial potential of mean torque V ( ψ ) = b sin ψ − k B T ln (cid:0) sin ψ (cid:1) . It is called the double-well potential in thepresent paper to oppose it to the ordinary Maier-Saupe one U MS ( ψ ) = b sin ψ .The barrier height is σ = b/ ( k B T ) − − ln [ b/ ( k B T )]. and the equation for the function η ( x ) (cid:16) − x (cid:17) [ η ′′ xx ( x ) − αxη ′ x ( x )] + 2 (cid:16) ρτ + 1 − α + 6 αx (cid:17) η ( x ) = 0 (12)We introduce a new variable y = x (13)The equation takes the form y ( y − η ′′ yy ( y ) + (cid:20) − αy + (cid:18) α + 12 (cid:19) y − (cid:21) η ′ y ( y ) − (cid:20) αy + 12 ( ρτ + 1 − α ) (cid:21) η ( y ) = 0 (14)It belongs to a class of the so-called confluent Heun’s equation [45]. Equation(14) has fundamental solutions that can be expressed via the confluent Heun’sfunction (CHF). The latter is a known special function [45], [47]. At presentit is realized explicitly in the only symbolic computational software packageMaple as HeunC and its derivative
HeunCP rime (see [47] for expert opin-ion on the merits and drawbacks of this computational tool in Maple). The8undamental solutions of (14) are η (1) ( y ) = HeunC (cid:18) − α, − / , − , − α , α − ρτ y (cid:19) (15) η (2) ( y ) = y / HeunC (cid:18) − α, / , − , − α , α − ρτ y (cid:19) (16)The spectrum of the eigenvalues for the parameter ρ is defined by the bound-ary conditions. For the latter we impose the usual requirement that the dis-tribution function f ( x, t ) must be finite at the boundary. In fact, due to thedivergence of our double-well Maier-Saupe uniaxial potential of mean torqueat the boundary ( V ( x ) → ∞ at x → ±
1) we require that f ( x, t ) must be zerothere lim x →± f ( x, t ) = 0 (17)The latter yields lim x →± η ( i ) ( x ) = 0 (18)for both i = 1 , lim y → η ( i ) ( y ) = 0 (19)We denote∆ n = ρ (1) n τ (20)Ω m = ρ (2) m τ (21)The equation for the spectrum of eigenvalues of ∆ n ( n = 1 , , , ... ) for η (1) n ( y )is HeunC (cid:18) − α, − / , − , − α , α − ∆ n y → (cid:19) = 0 (22)That for the spectrum of eigenvalues of Ω m ( m = 1 , , , ... ) for η (2) m ( y ) is HeunC (cid:18) − α, / , − , − α , α − Ω m y → (cid:19) = 0 (23)These equations can be solved only numerically but Maple easily copes withthis problem. One always obtains ∆ = 0 while the lowest Ω is always nonzero9 Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ (cid:144)H k B T L H D n L Ÿ D è è è è è è è è è è (cid:144)H k B T L H D n L è D ò ò ò ò ò ò ò ò ò ò (cid:144)H k B T L H D n L ò D ô ô ô ô ô ô ô ô ô ô (cid:144)H k B T L H D n L ô D ð ð ð ð ð ð ð ð ð ð (cid:144)H k B T L H D n L ð D Fig. 2. The dependence of the spectrum of the eigenvalues of ∆ n (obtained as thesolution of (22)) on the barrier height (∆ is the bottom line, ..., ∆ is the top one).This spectrum characterizes the first set of eigenfunctions (15) of the Smoluchowski’soperator L S if the Smoluchowski’s equation under study (3) to be written in theformal way ˙ f = L S f . and the corresponding term mainly determines the behavior of the distributionfunction f ( ψ, t ) and MFPT calculated with its help. The dependence of Ω n and ∆ m = 0 on the parameter α is depicted in Fig.2 and Fig.3 respectively.It should be stressed that at large barrier heights the value of Ω becomesvery small as Fig.3 testifies. We return to the variable ψ with the help of (4) and (13). Then the generalsolution of SE (3) can be written as f ( ψ, t ) = C η (1)1 ( ψ ) + ∞ X n =2 C n exp ( − ∆ n t/τ ) η (1) n ( ψ )+ ∞ X m =1 D m exp ( − Ω m t/τ ) η (2) m ( ψ ) (24)The crucial issue for the calculation of the coefficients C n and D m is thefollowing one: being the solution of the boundary problems both η (1) n ( ψ ) and10 Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ (cid:144)H k B T L - - H W m L Ÿ W ì ì ì ì ì ì ì ì ì ì (cid:144)H k B T L - - H W m L ì W è è è è è è è è è è (cid:144)H k B T L - - H W m L è W ò ò ò ò ò ò ò ò ò ò (cid:144)H k B T L - - H W m L ò W ô ô ô ô ô ô ô ô ô ô (cid:144)H k B T L - - H W m L ô W ð ð ð ð ð ð ð ð ð ð (cid:144)H k B T L - - H W m L ð W Fig. 3. The dependence of the spectrum of the eigenvalues of Ω m (obtained asthe solution of (23)) on the barrier height (Ω is the bottom line, ..., Ω is thetop one). This spectrum characterizes the second set of eigenfunctions (16) of theSmoluchowski’s operator L S if the Smoluchowski’s equation under study (3) to bewritten in the formal way ˙ f = L S f . η (2) m ( ψ ) are full sets of the orthogonal functions. Besides, as they belong todifferent fundamental solutions they are orthogonal to one another π Z dψ sin ψ η ( i ) l ( ψ ) η ( j ) m ( ψ ) ∝ δ lm δ ij (25)for i, j = 1 , η (1) n ( ψ ) or η (2) m ( ψ ). At time t = 0 (24) yields (with taking into account (2))2sin ψ δ ( ψ − ψ ) = ∞ X n =1 C n η (1) n ( ψ ) + ∞ X m =1 D m η (2) m ( ψ ) (26)Then the coefficients C n and D m are obtained by multiplying (26) by η (1) n ( ψ )or η (2) m ( ψ ) and integration over ψ with taking into account (25) C n = 2 η (1) n ( ψ ) π Z dψ sin ψ h η (1) n ( ψ ) i − (27)11nd D m = 2 η (2) m ( ψ ) π Z dψ sin ψ h η (2) m ( ψ ) i − (28)Thus, we have the explicit algorithm for obtaining the coefficients C n and D m along with the corresponding spectra ∆ n and Ω m that makes the solutionof our problem to be completed. Substitution of all these values into (24)yields the required probability distribution function. Of course, in practice atruncation of the series is necessary, i.e., replacement of the infinity by somefinite number M that is determined by the required accuracy. Below we provethat in practice the first terms M = 1 (i.e., those with the coefficients C and D ) provide sufficient accuracy (see Fig. 6 and discussion in Sec. 6). To exhibit how the result obtained may be useful we calculate with its helpthe dependence of MFPT on the barrier height α . Our effective Maier-Saupeuniaxial potential of mean torque is (see (7), (8) and (9))1 k B T V ( ψ ) = 2 α sin ψ − ln (cid:16) sin ψ (cid:17) (29)It is depicted in Fig. 1. It has two local minima the left of which is at ψ L = arcsin √ α ! (30)and the maximum at ψ M = π ψ = ψ L . It is worthy to recall thatour solution of SE f ( ψ, t ) is actually the conditional probability f ( ψ, t ) ≡ P ( ψ, t ; ψ ,
0) of finding the system at orientation ψ at time t , if the orientationwas ψ at time zero. Then, following Risken [6] we introduce the probabilityΩ( ψ , t ) of realizations which have started at ψ = ψ L and which have not yet12eached one of the boundaries ψ = 0 or ψ = ψ M up to the time t Ω( ψ , t ) = ψ M Z dψ sin ψ f ( ψ, t ) (32)The distribution function w ( ψ , T ) for the first passage time T (we use thenotations from [6] and note that this value is not to be confused with tempera-ture which will further enter only implicitly via α from (9) and via 1 / ( k B T ) =2 α/b ) is w ( ψ , T ) = − ∂ Ω( ψ , T ) ∂T (33)The moments of the first passage time distribution are T n ( ψ ) = ∞ Z dT T n w ( ψ , T ) (34)The value of interest for us here is the MFPT T T ( ψ ) = ∞ Z dT T w ( ψ , T ) (35)that characterizes the ability for the system to reach the barrier top ψ M start-ing from the left minimum ψ = ψ L . Thus, we have to calculate the quantity T ( ψ ) = − ψ M Z dψ sin ψ ∞ Z dT T ∂f ( ψ, T ) ∂T (36)After straightforward calculations we obtain1 τ T ( ψ ) = D Ω ψ M Z dψ sin ψ η (2)1 ( ψ )+ ∞ X n =1 ψ M Z dψ sin ψ (cid:20) C n ∆ n η (1) n ( ψ ) + D n Ω n η (2) n ( ψ ) (cid:21) (37)13 (cid:144)H k B T L H < T > (cid:144) Τ L - Kramers’estimite0 5 10 15 20 b (cid:144)H k B T L H < T > (cid:144) Τ L - Coffey et al. approach0 5 10 15 20 b (cid:144)H k B T L H < T > (cid:144) Τ L - present paper approach Fig. 4. The dependence of the mean first passage time for the transition from thewell of the symmetric double-well Maier-Saupe uniaxial potential of mean torque V ( ψ ) = b sin ψ − k B T ln (cid:0) sin ψ (cid:1) on the parameter b which characterizes the barrierheight. The dots are the result of our approach developed in the present paper. Thedashed line is the result of the Kramers’ estimate for the transition rate in the highfriction limit (45). The continuous line is the result of the approach developed byCoffey, Kalmykov, D´ejardin and their coauthors (42). Finally, we calculate the averaged MFPT < T > for the case when the initialcondition is thermodynamically averaged over the left well < T > = π Z dψ sin ψ exp − V ( ψ ) k B T ! − × ψ M Z dψ sin ψ exp − V ( ψ ) k B T ! T ( ψ ) (38)It is worthy to recall here that the coefficients C n and D n are the functionsof the initial condition ψ , i.e., explicitly C n ≡ C n ( ψ ) and D n ≡ D n ( ψ ).The results of numerical calculation of < T > are presented in Fig.4. Therethe results for the averaged MFPT < T > obtained from our numericalcalculations are compared with the analytical estimate from the Kramers’theory for the transition rate1 τ < T > ≈ K (39)14f we return again to the variable x = cos ψ and denote h ( x ) ≡ k B T V ( x ) = 2 α (1 − x ) − ln (1 − x ) (40)then the Kramers’ formula takes the form [19], [20],Γ K = q h ′′ xx ( x L ) | h ′′ xx ( x M ) | π exp {− [ h ( x M ) − h ( x L )] } (41)Also in Fig.4 the results of the approach of Coffey, Kalmykov, D´ejardin andtheir coauthors (CKD) are presented. The latter in our case of the potential(29) is eq. (1.18.1.7) of [26] < T >τ = 2 π/ Z dψ ′ e α sin ψ ′ − ln ( sin ψ ′ )sin ψ ′ ψ ′ Z dψ sin ψ e − α sin ψ + ln ( sin ψ )(42)Following CKD we introduce the longitudinal correlation function < cos ψ (0) cos ψ ( t ) >< cos ψ (0) > = π Z dψ sin ψ π Z dψ sin ψ cos ψ cos ψ f ( ψ, t ) × π Z dψ sin ψ π Z dψ sin ψ (cos ψ ) f ( ψ, t ) − (43)where the averaging < ... > is carried out with the help of the obtaineddistribution function f ( ψ, t ) given by (24). Also, we define the longitudinalcorrelation time τ || = ∞ Z dt < cos ψ (0) cos ψ ( t ) >< cos ψ (0) > (44)In Fig.5 the results of calculations of the value τ || are presented. They arecompared with the result of CKD that in the case of the potential (29) is eq.(2.10.25) of [26]1 τ τ CKD || = 2 Z < cos ψ > × π Z dψ ′ e α sin ψ ′ − ln ( sin ψ ′ )sin ψ ′ ψ ′ Z dψ sin ψ cos ψ e − α sin ψ + ln ( sin ψ ) (45)15 (cid:144)H k B T L H Τ þ (cid:144) Τ L - Coffey et al. approach0 5 10 15 20 b (cid:144)H k B T L H Τ þ (cid:144) Τ L - present paper approach Fig. 5. The dependence of the longitudinal correlation time for rotational mo-tion in the symmetric double-well Maier-Saupe uniaxial potential of mean torque V ( ψ ) = b sin ψ − k B T ln (cid:0) sin ψ (cid:1) on the parameter b which characterizes the bar-rier height. The dots are the result of our approach developed in the present paper.The line is the result of the approach developed by Coffey, Kalmykov, D´ejardin andtheir coauthors (45). where Z = π Z dψ sin ψ e − α sin ψ + ln ( sin ψ ) (46)and < cos ψ > = 1 Z π Z dψ sin ψ cos ψ e − α sin ψ + ln ( sin ψ ) (47)Finally, in Fig.6 we plot the dependence of our calculated MFPT on the trun-cation number M in the series (24). To clarify the terminology it should be stressed that we deal with barrierheights in the energetic units (i.e., conceive it to be the ratio of its value inthe natural units to the product of the Boltzmann constant and temperature).16
Ÿ Ÿ Ÿ Ÿ Ÿ H < T > (cid:144) Τ L Ÿ b (cid:144)H k B T L = è è è è è è H < T > (cid:144) Τ L è b (cid:144)H k B T L = ò ò ò ò ò ò H < T > (cid:144) Τ L ò b (cid:144)H k B T L = ô ô ô ô ô ô H < T > (cid:144) Τ L ô b (cid:144)H k B T L = ð ð ð ð ð ð H < T > (cid:144) Τ L ð b (cid:144)H k B T L = ì ì ì ì ì ì H < T > (cid:144) Τ L ì b (cid:144)H k B T L = Fig. 6. The dependence of the mean first passage time for the transition from thewell of the symmetric double-well Maier-Saupe uniaxial potential of mean torque V ( ψ ) = b sin ψ − k B T ln (cid:0) sin ψ (cid:1) on the truncation number M in the series (24)( b/ ( k B T ) = 1 is the bottom line, ..., b/ ( k B T ) = 20 is the top one). Thus, e.g., the low barrier regime may well be realized for high barriers in thenatural units at sufficiently high temperatures. Namely in this sense we furtheruse the notions of the intermediate to high barrier (low temperature) regionand the low barrier (high temperature) one.Fig.4 testifies our approach yields absolutely identical results with CKD one(42) in the intermediate to high barrier region, i.e., at barrier heights greaterthan b/ ( k B T ) ≈ σ ≈ b/ ( k B T ) greater than ≈
4. We conclude that in theintermediate to high barrier region our results are in full agreement with theknown literature ones.In the low barrier region our results yield noticeably smaller values for MFPTthan those predicted by both the Kramers’ estimate (41) and that of CKD(42). We conclude that in this region the formula (42) of CKD (or more ac-curately, eq. (1.18.1.7) of [26]) appreciably overestimates MFPT or in otherwords predicts too small values for the transition rates from the potentialwell. The reason of the discrepancy between our results for very small bar-rier heights for MFPT and those of CKD seen in Fig.4 is in the fact thattheir formula (42) (or more accurately, eq. (1.18.1.7) of [26]) is obtained inthe stationary limit. It seems that for very small barrier heights the tran-sient dynamics plays a crucial role and has to be taken into account explicitly.17he latter requirement is satisfied both in our approach and at deriving theCKD formula for the longitudinal correlation time (45) (or more accurately,eq. (2.10.25) of [26]). Fig. 5 testifies that for the longitudinal correlation timewe obtain absolute identity of our results with those of CKD (45) in the wholerange of barrier heights.The characteristic feature is that both MFPT (Fig.4) and the longitudinalcorrelation time (Fig. 5) are smaller than the relaxation time τ for very smallbarrier heights. It is noticeable that the formulas of CKD also lead to such aconclusion (see Fig.4 and Fig.5). The physical sense of it may be as follows. Therelaxation time τ is an average characteristic rather than a minimal possibletime of the processes in the system. Its value sets the time scale for the motionproceeding in the diffusional regime. The characteristic length figuring in therelaxation time is the radius of a rotating sphere that models the physicalobject of interest. Thus, in itself τ is a rather rough parameter that doesnot take into account either the characteristic distance of the motion or theform of the potential. When the barrier height in our double-well Maier-Saupepotential becomes very small the distance between its minimums also becomesvery small (they coalesce at ψ = π/
2) in contrast to the ordinary Maier-Saupeone where they still remain at ψ = 0 and ψ = π . In this case the distancefor the motion becomes commensurable or even less than the characteristiclength figuring in the relaxation time. In our opinion, it is not surprising that tomove a vanishingly small distance over a vanishingly small barrier the systemrequires lesser time than τ (MFPT < τ ).To prove the convergence of the obtained series (24) we plot in Fig.6 thedependence of the calculated value for MFPT on the truncation number M .This Fig. shows that M = 1 provides sufficient accuracy because the resultscease to change at further increase of M for all barrier heights. This fact can beanticipated for large barrier heights where the term M = 1 or more preciselythe term with the coefficient D overwhelmingly dominates the dynamics ofthe system due to very low values of the smallest non-vanishing eigenvalue Ω (see Fig. 3). However, our results testify that also for small barrier heights thethe terms with the coefficient C and D in the series (24) are sufficient.Our approach yields directly the probability distribution function rather thanthe statistical moments (i.e., values averaged over it). We recall that CKDenables one to obtain the reformulation of the problem as the differential-recurrence relations for the longitudinal and transverse correlation functions.A systematic way for their analysis has been suggested that yields ”exactanalytic solutions for the longitudinal and transverse complex susceptibilitiesand correlation times for many problems of practical interest” [26]. However,in some applications (e.g., to relaxation processes in NMR) it is desirable towork namely with the distribution function taken alone. The virtue of ourapproach is in providing such option.18ur approach yields the exact solution of SE (3) along with that of CKD. Letus formally write this equation as ˙ f = L S f where L S is the Smoluchowski’soperator. Then the difference of our exact solution from that of CKD can beexplained as follows. The authors of CKD expand f over Legendre polynomi-als while we expand it over the eigenfunctions of L S . The latter are expressedover CHF which at present can be considered as a well tabulated special func-tion. Efficient numerical algorithms for it are available in the literature letalone an explicit realization in Maple that makes its usage a routine prob-lem. In fact, to deal with CHF is no harder than with Legendre polynomials.Moreover, CKD approach leads to recurrence system for the coefficients inthe expansion when those with lower indexes are related to higher ones. Tocope with such system is a difficult problem in itself though it is successfullysolved in the above mentioned papers. In our exact solution the coefficientswith different indexes in the expansion are independent of each other. How-ever, besides the advantages of our solution over the approximate one there aredrawbacks. Within the framework of CKD approach the authors of the abovementioned papers manage to obtain an analytical expression for the smallestnon-vanishing eigenvalue. In our approach that is obtained as a solution of(23). The latter is to be solved numerically. In the present paper we do nottry to obtain an analytical expression for Ω from (23). Also, our approach isdeveloped only for symmetric uniaxial potential of mean torque. The crucialsubstitution (13) can work only for potentials with even powers of cos ψ orsin ψ and can not be generalized to odd powers leading to the asymmetry ofthe potential.In conclusion, we revisit a rather old problem of rotational motion in a double-well potential of mean torque. The corresponding physical process is withinthe framework of the escape rate theory that is more than seventy years old (ifto count from the cornerstone Kramers’ paper). Nevertheless, the some usefultools for its exact analytic treatment (the theory of CHF and its convenientnumerical realizations) were developed relatively recently. Their applicationprovides a new way to rederive old and well known results and yields a com-plementary development to them.Acknowledgements. The author is grateful to Dr. Yu.F. Zuev for helpful dis-cussions. The work was supported by the grant from RFBR. Maple is very efficient for numerical search of the roots of (22) and (23).Unfortunately, the calculations of the integrals containing
HeunC by Mapleproves to be too time consuming that makes the analysis of < T > withthe help of this program to be practically inconvenient. Besides, one can not19ontrol the processing at its usage [47]. For this reason, we prefer to employ theexplicit numerical realization of the CHF. We find that the required accuracycan well be retained at reducing the number of terms in the definition ofthe CHF by a series. As a consequence the processing time is diminishedsubstantially. For this purpose, we make use of the explicit representation ofthe CHF by a series [45]. It includes the three-term recurrence relation [45],[48] that in our case takes the form for the first solution of (14) η (1) ( y ) = HeunC (cid:18) − α, − / , − , − α , α − ρτ y (cid:19) = ∞ X k =0 v (1) k (cid:18) − α, − / , − , − α , α − ρτ (cid:19) y k (48)where A (1) k v (1) k = B (1) k v (1) k − + C (1) k v (1) k − (49)with the initial conditions v (1) − = 0, v (1)0 = 1. Here A (1) k = 1 − k (50) B (1) k = 1 + 4 α − k + 4 − α − ρτ k (51) C (1) k = − αk (cid:18) k − (cid:19) (52)For the second solution of (14) we have y − / η (2) ( y ) = HeunC (cid:18) − α, / , − , − α , α − ρτ y (cid:19) = ∞ X k =0 v (2) k (cid:18) − α, / , − , − α , α − ρτ (cid:19) y k (53)where A (2) k v (2) k = B (2) k v (2) k − + C (2) k v (2) k − (54)with the initial conditions v (2) − = 0, v (2)0 = 1. Here A (2) k = 1 + 12 k (55)20 (2) k = 1 + 4 α − k − ρτ k (56) C (2) k = − αk (57)To provide fast convergence of the series in (48) and (53) we truncate them,i.e., replace the infinity by some large but finite N which we choose from therequirement that the obtained values of < T > /τ cease to change withinthe limits of the required accuracy. Under these approximations we obtain anefficient numerical algorithm for calculating the integrals with the our CHF.21 eferences [1] N.V.Brilliantov, O.P. Revokatov, Molecular motion in disordered media,Moscow University Press, Moscow, 1996.[2] N.G. Van Kampen, Stochastic processes in physics and chemistry, 3-d ed.,Elsevier, 2007.[3] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry andthe natural sciences, 2-nd ed., Springer, 1985.[4] I. Oppenheim, K.E. Shuler, G.H. Weiss, Stochastic processes in chemicalphysics. The master equation. MIT Press, Cambrige, MA, 1977.[5] R.R. Vold, Deuterium NMR studies of molecular motion, in: Nuclear magneticresonance probes of molecular motion, ed. R. Tycko, Kluwer, 1994, 27 - 86.[6] H. Risken, The Fokker-Plank equation, 2-nd ed., Springer, 1989.[7] W.T. Coffey, D.A. Garanin, D.J. McCarthy, Crossover formulas in the Kramerstheory of thermally activated escape rates - application to spin systems, Adv.Chem. Phys. 117, Eds. I. Prigogine, S.A. Rice, John Wiley, 2001.[8] L.D. Favro, Theory of the rotational Brownian motion of a free rigid body,Phys.Rev. 119 (1960) 53-62.[9] J.H. Freed, Anisotropic rotational diffusion and electron spin resonancelinewidths, J.Chem.Phys. 41 (1964) 2077-2083.[10] P.L. Nordio, P. Busolin, Electron spin resonance line shapes in partially orientedsystems, J.Chem.Phys. 55 (1970) 5485-5490.[11] P.L. Nordio, G. Rigatti, U. Segre, Spin relaxation in nematic solvents,J.Chem.Phys. 56 (1972) 2117-2123.[12] C.F. Polnaszek, G.V. Bruno, J.H. Freed, ESR line shapes in the slow-motionalregion: Anisotropic liquids, J.Chem.Phys. 58 (1973) 3185-3199.[13] K.A. Valiev, E.N. Ivanov, Rotational Brownian motion, Uspechi Phys. Nauk(in Russian), 109 (1973) 31-64.[14] J.H. Freed, Stochastic-molecular theory of spin-relaxation for liquid crystals,J.Chem.Phys. 66 (1977) 4183-4199.[15] F.M. Kuni, A.A. Melikhov, B.A. Storonkin, Three dimensional rotationalrelaxation in external fields, Teor. Mat. Fiz., 34 (1978) 374-386.[16] C. Zannoni, A. Arcioni, P. Cavatorta, Fluorescence depolarization in liquidcrystals and membrane bilayers, Chem.Phys.Lipids 32 (1983) 179-250.[17] W.T. Coffey, M. Evens, P. Grigolini, Molecular diffusion and spectra, Wiley,1984.
18] R.R. Vold, R.L. Vold, Nuclear spin relaxation and molecular dynamics inordered systems: Models for molecular reorientation in thermotropic liquidcrystals, J.Chem.Phys. 88 (1987) 1443-1457.[19] H.A. Kramers, Brownian motion in a field of force and the diffusion model ofchemical reactions, Physica 7 (1940) 284-304.[20] P. H¨anggi, P. Talkner, M. Borkovec, Fifty years after Kramers’ equation:reaction rate theory, Rev.Mod.Phys. 62 (1990) 251-341.[21] R. Tarroni, C. Zannoni, On the rotational diffusion of asymmetric molecules inliquid crystals, J.Chem.Phys. 95 (1991) 4550-4564.[22] E. Berggren, F.L. Tarroni, and C. Zannoni, Rotational diffusion of uniaxialprobes in biaxial liquid crystal phases, J.Chem.Phys. 99 (1993) 6180-6200.[23] A. Polimeno, J.H. Freed, A many-body stochastic approach to rotational motionin liquids. In: Advances in chemical physics, ed. I. Prigogine and S.A. Rice, v.LXXXIII, Wiley, 1993.[24] W. Alexiewicz, Ensemble averages for Smoluchowski-Debye rotational diffusionin the presence of a two-angle-dependent reorienting force, Chem.Phys.Lett.320 (2000) 582-586.[25] B.U. Felderhof, Nonlinear response of a dipolar system with rotational diffusionto a rotating field, Phys.Rev. E66 (2002) 051503.[26] W.T. Coffey, Yu.P. Kalmykov, J.T. Waldron, The Langevin equation: withapplications to stochastic problems in physics, chemistry and electricalengineering, 2-nd ed., World Scientific, 2004.[27] W.T. Coffey, Y.P. Kalmykov, B. Ouari, S.V. Titov, Rotational diffusion andorientation relaxation of rodlike molecules in a biaxial liquid crystal phase,Physica A 368 (2006) 362-376.[28] Y.P. Kalmykov, S.V. Titov, W.T. Coffey, Inertial effects in the orientationalrelaxation of rodlike molecules in a uniaxial potential, J.Chem.Phys. 130 (2009)064110.[29] Y.P. Kalmykov, S.V. Titov, W.T. Coffey, Inertial and bias effects in therotational Brownian motion of rodlike molecules in a uniaxial potential,J.Chem.Phys. 134 (2011) 044530.[30] A.J. Martin, G. Meier, A. Saupe, Extended Debye theory for dielectricrelaxation in nematic liquid crystals, Symp. Faraday Soc., 5 (1971) 119-133.[31] B.A. Storonkin, Theory of dielectric relaxation in nematic liquid crystals,Kristallogr., 30 (1985) 841-853.[32] R.Y. Dong, Nuclear magnetic resonance of liquid crystals, Springer Verlag, 1997.[33] R.Y. Dong, Relaxation and the dynamics of molecules in the liquid crystallinephases, Progress in Nuclear Magnetic Resonance Spectroscopy 41 (2002) 115-151.
34] R.Y. Dong, Nuclear magnetic resonance spectroscopy of liquid crystals, WorldScientific Publishing Co., 2010.[35] W.F. Brown, Thermal fluctuations of a single-domain particle, Phys.Rev. 130(1963) 1677-1686.[36] W.F. Brown, Thermal fluctuations of fine ferromagnetic particles, IEEETransactions on Magnetics, 15 (1979) 1196-1208.[37] W.T. Coffey, Y.P. Kalmykov, Thermal fluctuations of magnetic nanoparticles:Fifty years after Brown, J. Appl. Phys. 112 (2012) 121301.[38] W.T. Coffey, D.S. F. Crothers, Yu.P. Kalmykov, E.S. Massawe, J.T. Waldron,Exact analytic formula for the correlation time of a single-domain ferromagneticparticle, Phys.Rev. E 49 (1994) 1869-1882.[39] W.T. Coffey, Finite integral representation of characteristic times oforientational relaxation processes: Application to the uniform bias force effectin relaxation in bistable potentials. Adv. Chem. Phys. 103 (1998) 259-333.[40] W.T. Coffey, D.S.F. Crothers, Comparison of methods for the calculation ofsuperparamagnetic relaxation times, Phys. Rev. E 54 (1996) 4768-4774.[41] W.T. Coffey, D.S.F. Crothers, S.V. Titov, Escape times for rigid Brownianrotators in a bistable potential from the time evolution of the Green functionand the characteristic time of the probability evolution, Physica A 298 (2001)330-350.[42] J. Xiang, J. Sun, N.S. Sampson, The importance of hinge sequence for loopfunction and catalytic activity in the reaction catalyzed by triosephosphateisomerase, J. Mol. Biol. 307 (2001) 1103-1112.[43] M. Kokkinidis, N.M. Glykos, V.E. Fadouloglou, Protein flexibility andenzymatic catalysis, Advances in protein chemistry and structural biology, 87(2012) 181-218.[44] R.W. Pastor, A. Szabo, Langevin dynamics of a linear rotor in a MaierSaupepotential: Kramers turnover of the flipping rate, J. Chem. Phys. 97 (1992) 5098-5112.[45] Heun Equations, A. Ronveaux (ed.), Oxford Univ. Press, 1995.[46] S.Y. Slavyanov, W. Lay, Special functions, a unified theory based onsingularities, Oxford: Oxford Mathematical Monographs, 2000.[47] P.P. Fiziev, D.R. Staicova, Solving systems of transendental equations involvingthe Heun functions, AJCM 2 (2012) 95-105.[48] P.P. Fiziev, Novel relations and new properties of confluent Heun functions andtheir derivatives of arbitrary order, J. Phys. A: Math. Theor. 43 (2010) 035203.34] R.Y. Dong, Nuclear magnetic resonance spectroscopy of liquid crystals, WorldScientific Publishing Co., 2010.[35] W.F. Brown, Thermal fluctuations of a single-domain particle, Phys.Rev. 130(1963) 1677-1686.[36] W.F. Brown, Thermal fluctuations of fine ferromagnetic particles, IEEETransactions on Magnetics, 15 (1979) 1196-1208.[37] W.T. Coffey, Y.P. Kalmykov, Thermal fluctuations of magnetic nanoparticles:Fifty years after Brown, J. Appl. Phys. 112 (2012) 121301.[38] W.T. Coffey, D.S. F. Crothers, Yu.P. Kalmykov, E.S. Massawe, J.T. Waldron,Exact analytic formula for the correlation time of a single-domain ferromagneticparticle, Phys.Rev. E 49 (1994) 1869-1882.[39] W.T. Coffey, Finite integral representation of characteristic times oforientational relaxation processes: Application to the uniform bias force effectin relaxation in bistable potentials. Adv. Chem. Phys. 103 (1998) 259-333.[40] W.T. Coffey, D.S.F. Crothers, Comparison of methods for the calculation ofsuperparamagnetic relaxation times, Phys. Rev. E 54 (1996) 4768-4774.[41] W.T. Coffey, D.S.F. Crothers, S.V. Titov, Escape times for rigid Brownianrotators in a bistable potential from the time evolution of the Green functionand the characteristic time of the probability evolution, Physica A 298 (2001)330-350.[42] J. Xiang, J. Sun, N.S. Sampson, The importance of hinge sequence for loopfunction and catalytic activity in the reaction catalyzed by triosephosphateisomerase, J. Mol. Biol. 307 (2001) 1103-1112.[43] M. Kokkinidis, N.M. Glykos, V.E. Fadouloglou, Protein flexibility andenzymatic catalysis, Advances in protein chemistry and structural biology, 87(2012) 181-218.[44] R.W. Pastor, A. Szabo, Langevin dynamics of a linear rotor in a MaierSaupepotential: Kramers turnover of the flipping rate, J. Chem. Phys. 97 (1992) 5098-5112.[45] Heun Equations, A. Ronveaux (ed.), Oxford Univ. Press, 1995.[46] S.Y. Slavyanov, W. Lay, Special functions, a unified theory based onsingularities, Oxford: Oxford Mathematical Monographs, 2000.[47] P.P. Fiziev, D.R. Staicova, Solving systems of transendental equations involvingthe Heun functions, AJCM 2 (2012) 95-105.[48] P.P. Fiziev, Novel relations and new properties of confluent Heun functions andtheir derivatives of arbitrary order, J. Phys. A: Math. Theor. 43 (2010) 035203.