Finite element solution of the Fokker-Planck equation for single domain particles
aa r X i v : . [ c ond - m a t . o t h e r] J a n Finite element solution of the Fokker-Planck equation for singledomain particles
N. V. Peskov
Faculty of Computational Mathematics and Cybernetics,Lomonosov Moscow State University, Moscow, Russian Federation. (Dated: January 7, 2020)
Abstract
The Fokker-Planck equation derived by Brown for the probability density function of the orienta-tion of the magnetic moment of single domain particles is one of the basic equations in the theoryof superparamagnetism. Usually this equation is solved by expanding the solution into a seriesof spherical harmonics, which in this case is a complex and cumbersome procedure. This articlepresents the implementation procedure and some results of the numerical solution of the Fokker-Planck equation using the finite element method. A method for creating a sequence of triangulargrids on the surface of a sphere based on an inscribed icosahedron is described. The equationsof the finite element method are derived and examples of numerical solutions are presented. Theprocesses of magnetization and demagnetization under heating of a particle with cubic magneticanisotropy are simulated. . INTRODUCTION A particle of ferromagnetic material below a certain critical size (typically 30 nm indiameter) constitutes a single domain particle meaning that it is in a state of uniformmagnetization for any applied field . The magnetic moment M per unit volume of such aparticle can be represented as a vector of constant magnitude, M = M s u , | u | = 1, where M s is the saturation magnetization per unit volume. In the course of time only the orientationof the magnetic moment, determined by the unit vector u , can change. The orientation isinfluenced by internal magnetocrystalline anisotropy, an external magnetic field, and randomfluctuations caused by thermal agitation.The magnetic properties of single domain particles in the absence of thermal fluctuationsdescribed by Stoner and Wohlfarth . The theory of the thermal fluctuations of the magne-tization of single domain particles was began with work of N´eel and was further developedby Brown .A particle with orientation u = (sin θ cos φ, sin θ sin φ, cos θ ) in a Cartesian coordinatesystem, where θ is the polar angle and φ is the azimuth, is assumed to be in internalthermodynamic equilibrium at temperature T , with Helmholtz free energy per unit volume E a ( θ, φ, T ). The particle is not necessarily in external equilibrium with the applied field H = H h , | h | = 1. The Gibbs free energy per unit volume is V ( θ, φ, T, H ) = E a ( θ, φ, T ) − M s H ( u · h ), which will be written below as V ( θ, φ ).In the absence of thermal agitation, changes of u are assumed to obey the Gilbertequation d u dt = − h ′ M s α ( u × ∇ V ) + h ′ ( u × ( u × ∇ V )) , (1)where t is the time, α is a dimensionless damping coefficient and h ′ = ( αγ ) / ((1 + α ) M s ), γ is the ratio of magnetic moment to angular momentum, and ∇ is the angular part of thegradient.In Ref. 4 the evolution of the magnetic moment was considered as a Brownian motionalong the surface of a unit sphere of a point, representing the orientation of the magneticmoment, subjected to the applied field and magnetic anisotropy. As the Langevin equa-tion for this motion Brown took Gilbert’s equation, supplemented by a random Gaussianwhite noise field, which takes into account the collisional damping. Using the obtainedLangevin equation, Brown derived the Fokker-Planck equation (FPE) for the probability2ensity function W ( θ, φ, t ) of orientations of magnetic moments, i.e. representing points onthe unit sphere.The FPE derived by Brown can be written in the form of a continuity equation ∂W∂t − ∇ · (cid:18) k ∇ W + d u dt W (cid:19) = 0 , (2)where the coefficient k should be chosen so that the Boltzmann distribution W B ∝ exp( − ( vV ) / ( k B T )) is a stationary solution of (2) for a particle of volume v , k B is theBoltzmann constant.Brown’s approach to the theory of magnetism of single domain particles opened up thepossibility of applying the methods developed in the theory of Brownian motion to thestudy of a superparamagnetism. A comprehensive review of these methods, as well as themost important results obtained with their help, are presented in the book by Coffey andKalmykov. Since the present paper is devoted to the numerical solution of Brown’s FPE, werestrict ourselves to a brief overview of commonly used methods for solving this equation.The solution of the equation (2) is usually associated with a decomposition of W in abasis consisting of angular eigenfunctions of the Laplace operator (spherical harmonics),which results in an infinite system of differential-recurrent equations for the coefficients ofdecomposition. The procedure for deriving of this system from the FPE equation is givenin Ref. 7. The system of differential-recurrent equations has the form ddt X ( t ) = A X ( t ) , (3)where X = { x ( t ) , x ( t ) , . . . } is the infinite vector of expansion coefficients, and A is amatrix that can depend on time.One of the peculiarities of the system (3) is that if the elements of the vector X are prop-erly ordered, then the matrix A becomes the d -diagonal matrix for any magnetic anisotropythat can be expressed as the finite combination of spherical harmonics. The number of di-agonals, d , is determined by the type of magnetic anisotropy. For isotropic particles, d = 3,for anisotropic particles, d > A may, in particular, arise due to the time dependenceof the applied field H . In studies related to the simulation of dynamic magnetic hysteresisor the calculation of dynamic magnetic susceptibilities a periodic applied field (ac field), H ( t ) = H cos ωt is usually considered. This field generates a time dependence through3os ωt of some elements of A . It seems that, for the first time, a study of dynamic hysteresisusing the numerical solution of FPE, transformed in the form (3), was undertaken in Ref. 8,where the hysteresis of isotropic superparamagnets was studied.For anisotropic particles for solving FPE under ac field the system (3) is reduced to alinear algebraic system for the coefficients F nm by substituting x m ( t ) = ∞ X n = −∞ F nm e ınωt , m = 0 , , . . . . The obtained linear system can be solved with the matrix sweep algorithm , or by the matrixcontinued fraction method .In connection with linear systems, continued fractions appear, in particular, when solv-ing a linear system with a 3-diagonal matrix by successively eliminating unknowns (Gaussmethod). For anisotropic particles, the equation system for the coefficients F nm has a d -diagonal matrix, d >
3. However, by grouping the unknowns into sub-vectors of the samelength so that for sub-vectors it is possible to obtain a linear system with a 3-diagonal ma-trix, whose elements will be matrices of small dimension . Therefore the solution of thenew system can be obtained as a matrix-valued continued fraction .The process of magnetic relaxation of single-domain particles can be described by theFPE with a constant applied field H = H (dc field). For a dc field, the system (3) can alsobe reduced to a linear algebraic system using the Laplace transform of X ( t ) and solved bythe matrix continued fraction method. Over the past two decades, many physical problemsassociated with single domain particles have been solved using the matrix continued fractionmethod. Various examples of such problems can be found in Ref. 6.An alternative method for investigating the statistical properties of single-domain parti-cles is the Monte Carlo method , which allows one to obtain macroscopic observables byaveraging microscopic ones. One of the difficulties in applying the Monte Carlo method tothe study of superparamagnetism is the uncertainty of the time scale of the Monte Carlosteps. Using a numerical FPE solution is expected to help overcome this difficulty .This work continues and develops the theme begun in the previous paper and devotedto the application of the finite element method (FEM) for solving Brown’s FPE. The FEMapproach to FPE is relatively simple and independent of the type of magnetic anisotropy.FEM directly gives the FPE solution, not spherical harmonics, as the matrix continuedfractions gives. The probability density function provided by FEM solution of FPE enables4ne to calculate any statistical characteristics of the single domain particle magnetization.This work uses the deterministic procedure for creating a triangular grid on the surfaceof a sphere based on an inscribed icosahedron, which is much simpler and more efficientthan the random number procedure described in the previous article . In Ref. 15, a FPEsolution was demonstrated with an applied ac field simulating dynamic magnetic hysteresis.Here are the solutions of the PFE with a dc field simulating the magnetization of a particleand its demagnetization with increasing temperature. The examples are calculated for cubicmagnetic anisotropy taking into account two anisotropy constants. II. FINITE ELEMENT SCHEME
Substituting d u /dt from Eq. (1) into Eq. (2), after some transformations, equation (2)can be written in a form convenient for applying the finite element method: ∂W∂τ = ∇ · (cid:20) ∇ W + W (cid:18) ∇ ˜ V + 1 α (cid:16) u × ∇ ˜ V (cid:17)(cid:19)(cid:21) , (4)where ˜ V = vV /k B T ; τ = t/ τ N is the dimensionless time, τ N = vM s (1 + α )2 k B T γα is the characteristic relaxation time.The next step of construction the FEM scheme is the generation of a triangular grid onthe surface of the sphere.
A. Triangular grid
In the present paper, for the finite element method, a regular triangular grid is constructedon the surface of the sphere. Here, regularity means that the positions of the grid nodesare calculated using a deterministic procedure and are not random, as was the case in theprevious work .An easy way to cover the sphere with triangles is to build a uniform triangular grid onthe surface of the inscribed icosahedron and transfer it to the surface of the sphere using thecentral projection. Such a grid can be called a ’raw grid’ . The faces of the icosahedron areequal equilateral triangles. Therefore, it is possible to build a uniform grid on the surface5f the icosahedron by dividing each of its edges into n equal segments and connecting thedividing points lying on adjacent edges with lines parallel to the closing edge. The result is auniform grid consisting of N t = 20 n triangles and N p = 10 n + 2 nodes. When one transfersthe grid to the surface of the sphere using the central projection, the initially uniform gridwill be distorted. The degree of grid distortion is usually characterized by two parameters:the ratio r of the lengths of the shortest and the longest linear elements, and the ratio r of the areas of the smallest and largest grid cells. For the ’raw grid’ with n = 81 (presentedbelow results were obtained with this n ) one has r = 0 . r = 0 . r and r . In particular, the authors of Ref. 16 report anoptimized grid with r = 0 .
786 and r = 0 .
952 for the number of nodes close to N p at n = 81. Since the main purpose of this paper is to demonstrate the potential of FEM toapply to the Brown equation, a very perfect grid that is difficult to build is not necessaryhere. Therefore, the grid used in the further calculations is constructed in the followingrather simple way. FIG. 1. Mapping a triangular grid from the icosahedron face to the surface of a sphere
The basic grid is the uniform grid on the surface of the icosahedron described above.Instead of the central projection, the following procedure is used to transfer the grid to thesphere. The vertices of the inscribed icosahedron remain in their places. The edges of theicosahedron, shown in Figure 1 by thick dashed lines, are mapped into arcs of large circlesconnecting adjacent vertices, which are shown by blue lines. Grid points on an arc, similarto grid points on an edge, divide each arc into n equal segments.6ach node of the base grid on the face of the icosahedron lies at the intersection of threestraight lines parallel to the edges that limit this face and passing through certain oppositenodes on the icosahedron edges. One of such nodes is shown in Fig. 1 by red dot. Theanalogues of these lines on the sphere are arcs of three great circles passing through relatednodes lying on arcs corresponding to edges. Arcs intersect in pairs, but all three do notintersect at one point. Therefore, the center of a spherical triangle with vertices located atthe points of pairwise intersection of the arcs is taken as the image on the surface of thesphere of the base grid node. This node is also indicated by red dot in Fig. 1. The gridconstructed in this way for n = 81 has the following distortion parameter values: r = 0 . r = 0 . FIG. 2. Triangular grid with n = 9, N p = 812, N t = 1620. Red dots indicate the vertices of theinscribed icosahedron. As an illustration of the grid constructed in this way, a triangular grid on a sphere at n = 9is shown in Fig. 2. Twelve nodes coinciding with the vertices of the inscribed icosahedron(marked in red in Fig. 2) each have 5 nearest neighbors. All other nodes have 6 nearestneighbors each. B. Finite element equations
The triangular grid constructed above can be considered as a polyhedron inscribed in thesphere with flat triangular faces and vertices P i , i = 1 , , . . . , N p , which are grid nodes. Let7 i be a neighborhood of the node P i , i.e. the union of triangular faces in which P i is one ofthe vertices. A neighborhood O i is composed of five adjacent triangles with a common vertex P i , if this vertex is the vertex of the inscribed icosahedron, or of six adjacent triangles, iftheir common vertex is not the vertex of the icosahedron. On the surface of the polyhedronwe define a real continuous function ϕ i , 0 ≤ ϕ i ≤
1, so that it is linear on each triangle, ϕ i ( P i ) = 1 and ϕ i ≡ O i . Functions ϕ i , i = 1 , , . . . , N p will be called finiteelements.Since there is a one-to-one mapping (central projection) between the polyhedron and thesphere, any function defined on the polyhedron can be considered on the sphere, and viceversa. Therefore, both sides of Eq. (4) can be multiplied by function ϕ i and integrated overthe surface of the sphere Z ϕ i ∂W∂τ d Ω = Z ϕ i ∇ · (cid:20) ∇ W + W (cid:18) ∇ ˜ V + 1 α (cid:16) u × ∇ ˜ V (cid:17)(cid:19)(cid:21) d Ω . Now we apply the Green formula to the integral on the right-hand side, in which there willbe no integral over the boundary, since the surface of the sphere has no boundary. Z ϕ i ∂W∂τ d Ω = − Z ∇ ϕ i · (cid:20) ∇ W + W (cid:18) ∇ ˜ V + 1 α (cid:16) u × ∇ ˜ V (cid:17)(cid:19)(cid:21) d Ω . (5)The solution to the last equation will be sought in the form W ( θ, φ, τ ) = N p X j =1 W j ( τ ) ϕ j ( θ, φ ) . (6)Substituting (6) into (5) one obtains the following system of linear ordinary differentialequations M ˙ W = − (cid:0) L + F + α − G (cid:1) W , (7)where W is the vector ( W , W , . . . , W N p ) T , dot denotes the derivative with respect to τ . M , L , F and G are square matrices of dimension N p × N p . Matrix elements are calculatedby the following formulas m ij = Z ϕ i ϕ j d Ω ,l ij = Z ( ∇ ϕ i · ∇ ϕ j ) d Ω , (8) f ij = Z ( ∇ ϕ i · ∇ ˜ V ) ϕ j d Ω ,g ij = Z ( ∇ ϕ i · ( u × ∇ ˜ V )) ϕ j d Ω .
8o calculate the matrix elements, the integrals over the sphere are approximated by theintegrals over the surface of the embedded polyhedron corresponding to the triangular grid.And by virtue of the definition of finite elements ϕ i , the integration domain is reduced to theintersection of neighborhoods O ij = O i ∩ O j . The set O ij at i = j consists of two adjacenttriangles with a common side P i P j provided that P i and P j are the nearest neighbors in thegrid.Matrices M and L depend only on the grid and are calculated accurately. m ij = S ( O ij )(1 + δ ij ) / , where S ( O ij ) is the area of O ij and δ ij is the Kronecker delta. The calculation of thegradient ∇ ϕ i is described in detail in Ref. 15. We only note here that since ϕ i is linear oneach triangle, therefore, its gradient is constant on each triangle, and the matrix elements l ij are also calculated accurately. Magnetic energy dependent matrix elements f ij and g ij can be calculated using numerical integration over grid triangles included in O ij . III. NUMERICAL EXAMPLES
To demonstrate the capabilities of the finite element method as applied to the Brownequation, we present results of two simulations. In both simulations, the parameter valuesfor Fe presented in Ref. 1 will be used. Fe possesses cubic magnetic anisotropy with internalmagnetic energy per unit volume that can by expressed as V a ( x, y, z ) = K [( xy ) + ( yz ) + ( zx ) ] + K ( xyz ) , (9)where K , K are the anisotropy constants and x, y, z, x + y + z = 1, are the guided cosinesof magnetic moment M , which can be considered as well as the Cartesian coordinates of apoint on the surface of a unit sphere. The function (9) defined on a sphere has 26 criticalpoints: 6 minima, 8 maxima and 12 saddle points. The directions corresponding to thesepoints will be used as the directions of the applied magnetic field in further simulation.For cubic anisotropy (9), the dimensionless function ˜ V can be written in spherical coor-dinates as ˜ V ( θ, φ )= ǫ a (cid:2) cot θ + (1 + κ cos θ ) cos φ sin φ (cid:3) sin θ − ǫ h [( h x cos φ + h y sin φ ) sin θ + h z cos θ ] , (10)9here κ = K /K , h x , h y , h z are the Cartesian components (guided cosines) of the unitvector h , and ǫ a = vK k B T , ǫ h = vM s Hk B T .
According to Ref. 1, K = 4.8 × erg/cm , K = 0.5 × erg/cm , therefore κ = 0 . M s at T = 20 ◦ C is equal 1714.0 emu/cm . Fordefiniteness, we consider a Fe particle of a cubic shape with an edge size 24 nm, therefore,the particle volume is v = 24 nm . Using the given parameter values, one can obtain ǫ a = 164 .
023 at T = 20 ◦ C (293 K). For simulation of magnetization the value of ǫ h = 4 ǫ a was taken.All calculations were performed on the grid described above for n = 81, the grid has N p = 65612 nodes and N t = 131220 triangles. Formula for the angular part of the gradienton a unit sphere ∇ = e θ ∂∂θ + e φ sin θ ∂∂φ , where e θ , e φ are the angular basic unit vectors of the spherical coordinate system, was usedto calculate ∇ ˜ V . Since component expression for u in spherical coordinates is u = (0 , , g ij were calculated by the formula g ij = Z O ij ∇ ϕ i · − θ ∂ ˜ V∂φ e θ + ∂ ˜ V∂θ e φ ! ϕ j d Ω . The integrals in f ij and g ij were estimated numerically for each triangle of the grid bydividing the triangle into 9 equal triangles and using the prismoidal formula.The first simulation presented here is connected with the process of magnetization ofa single-domain particle in a constant magnetic field. This problem was considered forvarious purposes in several papers , where FPE was solved using expansion in sphericalharmonics followed by the Laplace transform and the method of matrix continued fractions.We present only the results of the FEM FPE solution without any physical interpretation.FPE was solved for the initial equilibrium distribution W = Z − exp( − vV a /k B T ), Z = R exp( − vV a /k B T ) d Ω, for three directions of the applied field h corresponding to the criticalpoints V a . In Cartesian coordinates, these directions are expressed as h min = (0 , , h max = (1 / √ , / √ , / √
3) and h sad = (0 , / √ , / √ µ ( u · h ) minsadmax σ ( u · h ) τ δ W FIG. 3. Magnetization of a particle in constant field. (a) Average magnetization for three directionsof applied field. (b) Mean square root deviation of average magnetization. (c) Normalized deviationof current distribution W ( τ ) from final distribution W H . The value of damping parameter α = 0 . demonstration of the average magnetization µ ( u · h ), – the average projection of u on thedirection of applied field h : µ ( u · h ) = Z ( u · h ) W ( u .τ ) d Ω , the mean square root deviation of the magnetization: σ ( u · h ) = Z (( u · h ) − µ ( u · h )) W ( u , τ ) d Ω , and the normalized deviation of current density distribution from final equilibrium distribu-tion W H ∝ exp( − ˜ V ): δW = R | W ( τ ) − W H | d Ω R | W − W H | d Ω . All these integrals are calculated numerically using the numerical solution of FPE and theprismoidal formula. In particular, µ ( u · h ) ≈ N t X k =1 (( u k W k + u k W k + u k W k ) · h ) S k , u k , u k , u k are the radius-vectors of the vertices of k -th triangle, and S k is the areaof k -th triangle.Figure 3 shows three characteristics of magnetization process listed above obtained fromnumerical solution of FPE for initial equilibrium distribution at zero applied field. At turningon the external field the magnetic moment quickly alined with the direction of applied field.The mean square root of the average magnetization tends to zero, while the probabilitydensity distribution converges to the equilibrium (Boltzmann) distribution corresponding tothe applied field H .
200 400 600 800 1000 µ ( u · h ) minsadmax
200 400 600 800 1000 σ ( u · h ) o C200 400 600 800 1000 δ W
012 (a)(b)(c)
FIG. 4. Demagnetization of a particle under increasing temperature condition. The heating rate k T = 10 − K. (a) Average magnetization. (b) The standard deviation of the magnetization. (c)The deviation of current solution from the Boltzmann distribution at current temperature.
FPE was resolved for the damping parameter α of 0.01. It should be noted that the stiff-ness of the system of equations strongly depends on α and rapidly increases with decreasing α . The computation time also increases rapidly, since a very small time step is required toachieve the required accuracy. , The numerical examples presented here were calculated ona PC with a 4 GHz processor. With α = 1, the calculation takes less than 1 minute, with α = 0 . α = 0 . k T , T = T + k T τ .When ǫ h = 0 and the FPI numerical solution begins with the distribution W H – the finaldistribution in the first simulation, the distribution W ( τ ) changes very quickly and if α < α = 1 with initial condition W = W H and the heating rate k T = 10 − K per unit ofdimensionless time. The results, graphs of the average magnetization, the standard deviationof the magnetization, and the deviation of the current probability density distribution fromthe Boltzmann distribution at current T , in dependency on the particle temperature, areshown in Figure 4. Also, as in the first example, the graphs are plotted for the threedirections of the external field at which the particle was magnetized. Here the deviation δW is not normalized, δW ( τ ) = Z | W ( τ ) − W B ( T ( τ )) | d Ω . T, o C200 400 600 800 1000 µ ( u · h ) FIG. 5. Demagnetization of a particle under increasing temperature with different heating rates.Graphs of the average magnetization: blue line – k T = 10 − K, black line – k T = 10 − K, red line – k T = 10 − K. Solid lines – α = 1, dashed lines – α = 0 . When heated, the average magnetization of the magnetized particle remains almost at theinitial level until a certain temperature is reached, and then quickly drops. The temperatureat which the falling magnetization becomes less than a certain value depends on the heatingrate. For example, the average magnetization of a particle magnetized in an external fielddirected to the saddle point of V a , becomes less than 0.01 at T = 940 ◦ C if k t = 10 − K, at13 = 580 ◦ C if k t = 10 − K, and at T = 390 ◦ C if k t = 10 − K. These values are obtainedfrom the FPE solution with α = 1. With a decrease in α , these temperatures are alsodecrease. In particular, for α = 0 . T = 820K, T = 510K, and T = 340K. As mentionedabove, solving Eq. (7) with the initial condition W H for α = 0 . α = 1 at thetime τ = 1. The graphs of the average magnetization versus temperature for three differentheating rates are shown in Figure 5. IV. CONCLUSION
The results of this and previous works show that the finite element method allowsone to efficiently solve Brown’s FPI for single-domain particles with apparently much lowercomputational and programmatic efforts than the traditional method of decomposition intospherical harmonics.The procedure for generating triangular grids on the surface of a sphere using an inscribedicosahedron enables one to create an infinite sequence of regular (deterministic) grids with afairly high quality. A comparison of the results obtained on grids with an increasing numberof nodes makes it possible to draw a conclusion about the convergence of the numericalmethod. In this work, the calculations were performed on grids with n = 72, 81, and99. The results of calculations on these grids almost coincided, which makes it possibleto conclude that it is inappropriate to increase the number of nodes to solve the problemsconsidered.All matrices of the system (7) have the same structure, nonzero elements in all matricesare in the same places. The number of nonzero elements, N nz = 7 N p −
12, is relatively small,so storing matrices does not require a lot of memory, and many calculations can be doneon a regular PC. Another feature of the system (7) is that it is not resolved with respect totime derivatives. Therefore, for its numerical solution, it is necessary to use special codesfor linear implicit systems, for example, the LSODIS codes. All the numerical examplespresented in the paper were run on a regular PC using the MATLAB environment. Tosolve the system (7), the MATLAB function ’ode15s’ was used. However, as parameter α decreases, the system stiffness increases rapidly, and solving a problem on a PC becomesproblematic due to a very long computational time. Therefore, it is advisable to develop14pecial codes for solving the system (7), taking into account its features and capabilities ofmulti-core processors. B.D. Cullity, C.D. Graham, Introduction to magnetic materials, IEEE Press, 2009. E.C. Stoner, E.P. Wohlfarth, A mechanism of magnetic hysteresis in heterogeneous alloys,Philos. Trans. R. Soc. London, Ser. A, (1948) 599-642. L. N´eel, Influence des fluctuations thermiques sur laimantation de grains ferromagn´etiques trsfins, C. R. Acad. Sci. Paris, (1949) 664-666. W.F. Brown Jr., Thermal fluctuations of a single-domain particle, Phys. Rev., (1963)1677-1686. T.L. Gilbert, A phenomenological theory of damping in ferromagnetic materials, IEEE Trans.Magn., (2004) 3443-3449. DOI: 10.1109/TMAG.2004.836740. W.T. Coffey, Yu.P. Kalmykov, The Langevin Equation, 4nd ed. World Scientific, Singapore,2017. Yu.P. Kalmykov, S.V. Titov, Derivation of matrix elements for the system of moment equationsgoverning the kinetics of superparamagnetic particles, J. Magn. Magn. Mater., (2000) 233-243. DOI: 10.1016/S0304-8853(99)00594-6 V.A. Ignatchenko, B.S. Gekht, Dynamic hysteresis of a superparamagnet, Zh. Eksp. Teor. Fiz.,67 (1974) 1506-1515. [Sov. Phys. JETP, (1975) 750-758. I.S. Poperechny, Yu.L. Raikher, V.I. Stepanov, Dynamic magnetic hysteresis in single-domainparticles with uniaxial anisotropy, Phys. Rev. B, (2010) 174423. DOI: 10.1103/Phys-RevB.82.174423 H. Risken, The FokkerPlanck Equation, Methods of Solutions and Applications, Springer-Verlag, Berlin, 1989. H. Zhaoa, G. Zhu, Matrix-valued continued fractions, J. Approx. Theory, (2003) 136-152.DOI: 10.1016/S0021-9045(02)00016-3 D.A. Dimitrov, G.M. Wysin, Magnetic properties of superparamagnetic particles by a MonteCarlo method. Phys. Rev. B, (1996) 9237-9241. DOI: 10.1103/PhysRevB.54.9237 X.Z. Cheng, M.B.A. Jalil, H.K. Lee, Y. Okabe, Mapping the Monte Carlo Scheme to LangevinDynamics: A Fokker-Planck Approach. Phys. Rev. Lett., (2006) 067208. DOI: 10.1103/Phys- evLett.96.067208 P. V. Melenev, Yu. L. Raikher, V. V. Rusakov, R. Perzynski, Time quantification for MonteCarlo modeling of superparamagnetic relaxation. Phys. Rev. B (2012) 104423. DOI:10.1103/PhysRevB.86.104423 N.V. Peskov, N.L. Semendyaeva, Numerical solution of Fokker-Planck equation forsingle domain particles. Phusica B: Condensed Matter, (2019) 142148. DOI:10.1016/j.physb.2019.07.004 R.P. Heikes, D.A. Randall, C.S. Konor, Optimized Icosahedral Grids: Performance of Finite-Difference Operators and Multigrid Solver. Mon. Weather Rev., (2013) 44504469. DOI:10.1175/MWR-D-12-00236.1 W.T. Coffey, D.S.F. Crothers, Yu.P. Kalmykov, J.T. Waldron, Constant-magnetic-field effect inN´eel relaxation of single-domain ferromagnetic particles, Phys. Rev. B, (1995) 15947-15956.DOI: 10.1103/PhysRevB.51.15947 Yu.P. Kalmykov, S.V. Titov, W.T. Coffey, Longitudinal complex magnetic susceptibility andrelaxation time of superparamagnetic particles with cubic magnetic anisotropy, Phys. Rev. B, (1998) 3267-3276. DOI: 10.1103/PhysRevB.58.3267 W.T. Coffeya, Yu.P. Kalmykovb, S.V. Titov. Nonlinear response of fine superparamagneticparticles to the sudden change of a strong uniform DC magnetic field, J. Magn. Magn. Mater., (2002) 400414. DOI: 10.1016/S0304-8853(01)00951-9 Yu.P. Kalmykov, W.T. Coffey, B. Ouaria, S.V. Titov, Damping dependence of the magnetizationrelaxation time of single-domain ferromagnetic particles, J. Magn. Magn. Mater., (2005)372384. DOI: 10.1016/j.jmmm.2004.11.233 B. Ouari, Yu.P. Kalmykov, Dynamics of the magnetization of single domain particles havingtriaxial anisotropy subjected to a uniform dc magnetic field, J. Appl. Phys., (2006) 123912.DOI: 10.1063/1.2399304 B. Ouari, S. Aktaou, Yu.P. Kalmykov, Reversal time of the magnetization of antiferromagneticnanoparticles, Phys. Rev. B, (2010) 024412. DOI: 10.1103/PhysRevB.81.024412(2010) 024412. DOI: 10.1103/PhysRevB.81.024412