format: text

text:
# Linear and nonlinear superparamagnetic relaxation at high anisotropy barriers

## I Introduction

The problem of superparamagnetic relaxation in single-domain ferroparticles formulated, explained and basically analyzed by NéelNeel-AGCR:49  about fifty years ago, since then keeps to attract attention. Nowadays this interest is mainly due to the expanding number of nanometer granular magnetic media used in information storage and related high technologies.

When analyzing magnetic dispersions, solid or fluid, a promising idea is to evaluate the granulometric content, particle material parameters and relaxation rates by combining the data on linear and nonlinear dynamic susceptibilities. Since recently, this approach (it originates from the spin glass science) became quite feasible in experimental realization.BiOh-JPSJ:96  However, to benefit from it, one needs an adequate model. Surprisingly, until nowadays the NéelNeel-AGCR:49  concept of superparamagnetic behavior of fine magnetic particles that had been substantially advanced by BrownBrown-JAP:59 ; Brown-PR:63  and refined by numerous researchers (see the review articleDoFi-ACP:97  with about 400 references), lacks a nonlinear extension.

In Ref. RaSt-PRB:97,  we begun to fill up this gap and proposed a numerical procedure involving continuous fractions by means of which the linear and cubic susceptibilities for a solid system of uniaxial fine particles could be obtained. With allowance for the polydispersity of real samples, the worked out description provided a fairly good agreement with the dynamic magnetic measurements taken on Co-Cu nanocomposites.BiOh-JPSJ:96  Recently, our approach was used successfullySpFi-JMMM:01  for the linear and cubic susceptibilities of the samples of randomly oriented $`\gamma `$-Fe2O3 nanoparticles. Hereby we carry on the build-up of the nonlinear superparamagnetic relaxation theory by working out a set of compact and accurate analytical expressions that considerably facilitate calculations as well as experiment interpretation.

The paper is arranged in the following way. In Sec. II we discuss the problem of superparamagnetic relaxation and show the way to obtain the asymptotic solution for the micromagnetic Fokker-Planck equation in the uniaxial case. In Sec. III the perturbative expansions for the orientational distribution function are obtained, which are used in Sec. IV to construct asymptotic expressions for the nonlinear dynamic susceptibilities. The explicit forms of those expansions are given and their accuracy is proved by comparison with the results of numerical calculations. Sec. V contains the enveloping discussion.

## II Superparamagnetic relaxation times

### II.1 Uniaxially anisotropic particle

The cornerstone of the superparamagnetic relaxation theory is the Arrhenius-like law for the relaxation rate of a magnetic moment of a single domain particle predicted by Néel in 1949. The framework of this classical problem is as follows. Consider an immobile (e.g. fixed inside a solid matrix) single-domain grain of a volume $`v`$. This particle possesses a uniaxial volume magnetic anisotropy, $`K`$ being its energy density and $`𝒏`$ its easy axis direction. Since the temperature $`T`$ is assumed to be much lower than the Curie point, the particle magnetization $`I`$, as a specific parameter, is practically constant and the magnitude of the particle magnetic moment may be written as $`\mu =Iv`$. Denoting its direction by a unit vector $`𝒆`$, one concludes that the magnetic state of such a particle is exhaustively characterized by a pair of vectors: $`𝝁=Iv𝒆`$ and $`𝒏`$. Thence, the orientation-dependent part of the particle energy (in the absence of external magnetic fields) is

$$U=Kv(𝒆𝒏)^2,$$
(1)
where $`K`$ is assumed to be positive. As Eq. (1) shows, this energy has two equal minima. They are separated by the potential barrier of the height $`Kv`$ and correspond to $`𝒆\pm 𝒏`$ because for the magnetic moment $`𝒆`$ the directions $`𝒏`$ and $`\mathbf{}𝒏`$ are equivalent. At zero temperature, the magnetic moment $`𝒆`$, once located in a particular potential well, is confined there forever. At finite temperature, the probability of an overbarrier (interwell) transition becomes non-zero. If the ratio $`\sigma Kv/kT`$ is high enough, the transition rate is exponential thus yielding the Néel law $`\tau \mathrm{exp}(\sigma )`$ for the reference time $`\tau `$ of the particle re-magnetization.

BrownBrown-PR:63  shaped up those semi-qualitative considerations into a rigorous Sturm-Liouville eigenvalue problem by deriving the micromagnetic kinetic equation

$$2\tau _DW/t=\widehat{𝑱}W\widehat{𝑱}(U/kT+\mathrm{ln}W),$$
(2)
where $`W(𝒆,t)`$ is the orientational distribution function of the magnetic moment, $`\widehat{𝑱}=(𝒆\times /𝒆)`$ is the infinitesimal rotation operator with respect to $`𝒆`$, and the time $`\tau _D`$ is introduced below by formula (4). Generally speaking, equation (2) is incomplete since a gyromagnetic term is absent there. This means that the consideration is limited by the frequency range $`\omega \tau _01`$, where $`\tau _0`$ is the relaxation time of the Larmor precession of the particle magnetic moment in the internal anisotropy field $`H_a=2K/I`$, where $`K`$ includes the possible shape contribution. Comparing this condition with the other one, $`\omega _L\tau _01`$, that evidences a low-to-moderate quality factor of the Larmor precession for real nanodisperse ferrites, one estimates the allowed frequency as $`\omega \omega _L`$ that means, in fact, a fairly wide range.111We remark that, in principle, there might occur a situation, where $`\tau `$ is very long, i.e., the quality factor of precession is very high. Then the gyromagnetic term in the kinetic equation must be retained so that Larmor precession begins to interact with the superparamagnetic (longitudinal) relaxation. An example of such a situation is considered in Ref. GaSv-PRL:00, , where some interesting nonlinear effects are found.

In the statistical description delivered by Eq. (2), the observed (macroscopic) magnetic moment per particle is given by the average

$$𝒎(t)=\mu 𝒆=𝒆W(𝒆,t)𝑑𝒆,$$
(3)
Note that with allowance for Eq. (1) the function $`W`$ has a parametric dependency on the vector $`𝒏`$ so that, in fact, the angular argument of $`W`$ is $`(𝒆𝒏)`$.

The magnetodynamic equation underlying the Brown kinetic equation (2) can be either that by Landau & Lifshitz or that by Gilbert. To be specific, we adopt the former one. Thence, the reference relaxation time in Eq. (2) writes as

$$\tau _D=Iv/2\alpha \gamma kT,$$
(4)
where $`\gamma `$ is the gyromagnetic ratio for electrons and $`\alpha `$ is the precession damping (spin-lattice relaxation) phenomenological parameter.

Assuming uniaxial symmetry of the time-dependent solution and separating the variables in Eq. (2) in the form

$$W(𝒆,t)=\frac{1}{2\pi }_{\mathrm{}=0}^{\mathrm{}}A_{\mathrm{}}\psi _{\mathrm{}}(𝒆𝒏)\mathrm{exp}(\lambda _{\mathrm{}}t/2\tau _D),$$
(5)
where the amplitudes $`A_{\mathrm{}}`$ depend on the initial perturbation, one arrives at the spectral problem

$$\widehat{L}\psi _{\mathrm{}}=\lambda _{\mathrm{}}\psi _{\mathrm{}},\widehat{L}\widehat{𝑱}\left[\mathrm{\hspace{0.17em}2}\sigma (𝒆𝒏)(𝒆\times 𝒏)\widehat{𝑱}\right],$$
(6)
where the non-negativity of the decrements $`\lambda _{\mathrm{}}`$ can be proven easily. Expanding the eigenmodes $`\psi _{\mathrm{}}`$ in the Legendre polynomial series

$$\psi _{\mathrm{}}=\frac{1}{2}\underset{k=1}{\overset{\mathrm{}}{}}(2k+1)b_k^{(\mathrm{})}P_k(\mathrm{cos}\theta ),k=\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}5},\mathrm{},$$
(7)
where $`\theta `$ is the angle between $`𝒆`$ and $`𝒏`$, one arrives at the homogeneous tridiagonal recurrence relation

$$\left[1\frac{\lambda _{\mathrm{}}}{k(k+1)}\right]b_k^{(\mathrm{})}2\sigma \left[\frac{k1}{(2k1)(2k+1)}b_{k2}^{(\mathrm{})}+\frac{1}{(2k1)(2k+3)}b_k^{(\mathrm{})}\frac{k+2}{(2k+1)(2k+3)}b_{k+2}^{(\mathrm{})}\right]=0,$$
(8)
Note that Eqs. (5)–(8) describe only the longitudinal (with respect to the easy axis) relaxation of the magnetic moment. We remark that under condition $`\omega \omega _\mathrm{L}`$, i.e., far from the ferromagnetic resonance range, the transversal components of $`𝒎=\mu 𝒆`$ are of minor importance.

### II.2 Interwell mode

Spectral equation (6) describes the temperature-induced (fluctuation) motions of the vector $`𝒆`$ in the orientational potential with a symmetrical profile (1). With respect to the time dependence, the set of possible eigenmodes splits into two categories: interwell (overbarrier) transitions and intrawell wanderings. In the spectral problem (6) the interwell transitions of the magnetic moment are associated with the single eigenvalue $`\lambda _1`$. As the rigorous analysis shows,Stor-K:85  it drastically differs from the others: whereas for $`\mathrm{}1`$ all the $`\lambda _{\mathrm{}}`$ gradually grow with $`\sigma `$, the decrement $`\lambda _1`$ exponentially falls down proportionally to $`\mathrm{exp}(\sigma )`$.

In the opposite limit $`\sigma 0`$, all the decrements, including $`\lambda _1`$, tend to the sequence $`\lambda _{\mathrm{}}=\mathrm{}(\mathrm{}+1)`$ and thus become of the same order of magnitude. This regime corresponds to a vanishing anisotropy so that the difference between the inter- and intrawell motions disappear, and the magnetic moment diffuses almost freely over all the $`4\pi `$ radians with the reference time $`\tau _D`$ introduced by Eq. (4).

From Eqs. (3) and (5) one finds that the longitudinal component of the magnetic moment evolves according to

$$m(t)=\mu \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}A_{\mathrm{}}e^{\lambda _{\mathrm{}}t/2\tau _D}_1^1x\psi _{\mathrm{}}𝑑x,$$
(9)
where $`x=\mathrm{cos}\theta =(𝒆𝒏)`$. For a symmetrical potential like (1) the equilibrium value $`m_0`$ of the particle magnetic moment is zero.

With the above-mentioned structure of the eigenvalue spectrum, the term with $`\mathrm{}=1`$ in Eq. (9), being proportional to $`\mathrm{exp}(e^\sigma t/\tau _D)`$, at $`\sigma >1`$ is far more long-living than any other one. The dominating rôle of the decrement $`\lambda _1`$ had been proven by Brown, and for it he had derivedBrown-PR:63  the asymptotic expression

$$\lambda _\text{B}=(4/\sqrt{\pi })\sigma ^{3/2}e^\sigma ,(\sigma 1).$$
(10)
In a short time after, using a continued fraction method, Aharoni constructedAhar-PR:64  for $`\lambda _1`$ a fairly long power series in $`\sigma `$ and also showed numerically that Brown’s expression (10) resembles the exact one with the accuracy of several percent for $`\sigma 3`$. In the 90’s the eigenvalue $`\lambda _1`$ became a subject of extensive studies. Efficient numerical procedures were developed Kalm-PRE:00  and a number of extrapolation formulas with a good overall accuracy were proposedBeBe-PRB:92 ; Ahar-PRB:92 ; CoCr-JMMM:94 ; CrCr-JAP:94 .

### II.3 Asymptotic solution of the Brown equation

The study that we describe below was inspired by our work on fitting the dynamic susceptibilities measurements for real assemblies of fine particles. Those data typically describe polydisperse systems in the low-frequency bandwidth $`\omega /2\pi =1`$–103 Hz. As $`\tau _010^9`$ s or smaller, then, using formula (10) for estimations, one concludes that the mentioned frequency interval becomes a dispersion range for the interwell (superparamagnetic) mode at

$$\omega \tau _0e^\sigma 1\mathrm{that}\mathrm{is}\sigma 10.$$

For temperatures up to 300 K this condition holds for quite a number of nanomagnetic systems.

Application of the best fit procedure to a set of experimental data implies numerous re-calculations of the linear and non-linear susceptibility curves $`\chi ^{(k)}`$ of the assembly. Any such curve, due to a considerable polydispersity of the particles, is a superposition of a great number of partial curves $`\chi ^{(k)}(\sigma )`$ spread over a wide size (or, in the dimensionless form, $`\sigma `$) range. For successful processing, one needs a fast and very accurate algorithm to evaluate $`\chi ^{(k)}(\sigma )`$ everywhere including the domain $`\sigma 1`$. The existing extrapolation formulas are no good for that purpose due to their ill-controllable error accumulation. A plausible way out is an asymptotic in $`\sigma ^1`$ solution of Eq. (6). In the course of the fitting procedure, this approximation can be easily matched in the intermediate $`\sigma `$ range with the well-known expansions for the small $`\sigma `$ end.

It is noteworthy that some 20 years ago Brown himself resumedBrown-P:77 ; Brown-IEEE:79  studies on $`\lambda _1`$ and modified the pre-exponential factor in Eq. (10) transforming it into an asymptotic series in $`\sigma ^1`$. On the base of Eq. (6) he had constructed an integral recurrence procedure, and evaluated $`\lambda _1`$ down to terms $`1/\sigma ^{10}`$. What we do below, is, in fact, carrying on this line of analysis that had not been touched since then. Our method advances Brown’s results in two aspects. First, for $`\lambda _1`$ it is more simple. Second, it provides not only the eigenvalue but the eigenfunction as well. Only having the latter in possession, one is able to obtain theoretical expressions for the directly measurable quantities that is the susceptibilities $`\chi ^{(k)}`$.

Taking Eq. (6) as the starting point, we remark its equilibrium solution

$`\psi _0=Z_0^1\mathrm{exp}(\sigma x^2),Z_0=2R(\sigma ),`$ (11)
$`R(\sigma )={\displaystyle _0^1}\mathrm{exp}(\sigma x^2)𝑑x,`$
that corresponds to $`\mathrm{}=0`$ and $`\lambda _0=0`$ and note the asymptotic expansion for the partition integral $`R(\sigma )`$ found in Ref. RaSh-JETP:74, :

$$R(\sigma )=e^\sigma G/2\sigma ,G(\sigma )1+\frac{1}{2\sigma }+\frac{3}{4\sigma ^2}+\frac{15}{8\sigma ^3}+\mathrm{}+\frac{(2n1)!!}{2^n\sigma ^n}+\mathrm{}.$$
(12)
The operator $`\widehat{L}`$ in Eq. (6) is not self-conjugated and thus produces two sets of eigenfunctions, which obey the respective equations

$$\widehat{L}\psi _k=\lambda _k\psi _k,\widehat{L}^+\phi _j=\lambda _j\phi _j;$$
(13)
here $`+`$ denotes Hermitian conjugation. The eigenfunctions of these two families are orthonormalized and related to each other in a simple way:

$$\psi _k=\psi _0\phi _k,_1^1𝑑x\phi _j\psi _k=\delta _{jk}.$$
(14)
Qualitatively, from Eq. (14) one may say that $`\phi _k`$ are the same eigenfunctions but “stripped” of the exponential equilibrium solution $`\psi _0`$. Substituting Eq. (14) in Eq. (6), one gets two useful relationships

$$\widehat{𝑱}\psi _0\widehat{𝑱}\phi _k=\lambda _k\psi _0\phi _k,\psi _0(\widehat{𝑱}\phi _j)(\widehat{𝑱}\phi _k)𝑑x=\lambda _k\delta _{jk},$$
(15)
where the second one follows from the first after multiplication by $`\phi _i`$ and integration by parts. Note that in the second formula action of each operator reaches no farther than the nearest closing parenthesis.

On rewriting Eq. (15.1) in terms of a single orientational variable $`x=(𝒆𝒏)`$, the spectral problem takes the form

$$\frac{d}{dx}\left[\psi _0(1x^2)\frac{d\phi _k}{dx}\right]=\lambda _k\psi _0\phi _k.$$
(16)
In the equilibrium state Eq. (16) reduces to

$$\frac{d}{dx}\left[\psi _0(1x^2)\frac{d\phi _0}{dx}\right]=0,$$
(17)
whose normalized solution is $`\phi _0=1`$. This solution, being a true equilibrium one, turns the inner part of the brackets, i.e., the probability flux in the kinetic equation (2), into identical zero.

As remarked in Sec. II.2, at $`\sigma 1`$ the most long-living non-stationary solution of Eq. (16) is the eigenfunction with $`\mathrm{}=1`$, whose eigenvalue is exponentially small, see Brown’s estimation (10). We use this circumstance for approximate evaluation of $`\phi _1`$ in the $`\sigma 1`$ limit by neglecting the right-hand side of Eq. (16) for $`\mathrm{}=1`$. On doing that, the equation obtained for the function $`\phi _1`$ formally coincides with equation (17) for $`\phi _0`$. However, the essential difference is that now the content of the bracket is non-zero:

$$\psi _0(1x^2)\frac{d\phi _1}{dx}=\frac{1}{2}C,$$
(18)
where $`\frac{1}{2}C`$ is the integration constant. Note also that, contrary to $`\phi _0`$, the sought for solution $`\phi _1`$ is odd in $`x`$.

Using the explicit form of $`\psi _0`$ from Eq. (11) and integrating, one gets for $`x>0`$

$$\phi _1=CR_0^x\frac{e^{\sigma x^2}}{1x^2}𝑑x=CR_0^xe^{\sigma x^2}(1+x^2+x^4+x^6+\mathrm{})𝑑x.$$
(19)
The integrals in expansion (19) are akin. Denoting

$$F_n=_0^xx^{2n}e^{\sigma x^2}𝑑x,$$

one can easily write for them the recurrence relation and “initial” condition as

$$F_n=\frac{}{\sigma }F_{n1},F_0=\frac{\sqrt{\pi }}{2\sqrt{\sigma }}\mathrm{erf}(\sqrt{\sigma }x),$$
(20)
respectively. Using the asymptotics of the error integral, with the exponential accuracy in $`\sigma `$ one finds

$$F_n=[(2n1)!!/2^n\sigma ^n]F_0,F_0\sqrt{\pi }/2\sqrt{\sigma }.$$
(21)
Comparing this with expression (12) for the function $`G`$, we get the representation

$$\phi _1(x>0)CRF_0G.$$
(22)
Applying to Eq. (22) the normalizing condition (14), one evaluates the constant as $`C=1/RF_0G`$. Therefore, from Eqs. (20)–(22) the principal relaxational eigenmode determined with the $`\mathrm{exp}(\sigma )`$ accuracy emerges as an odd step function

$$\phi _1(x)\{\begin{array}{cc}1\hfill & \text{for}x<0,\hfill \\ 1\hfill & \text{for}x>0.\hfill \end{array}$$
(23)
In Fig. 1 the limiting contour (23) is shown against the exact curves $`\phi _1(x)`$ obtained by solving numerically Eq. (8) for several values of $`\sigma `$. We remark that in the statistical calculations carried out below the typical integrals are of two kinds. In the first, the integrand consists of the product of $`\phi _1\psi _0`$ and some non-exponential function. As $`\psi _0\mathrm{exp}\sigma x^2`$, the details of behavior of $`\phi _1`$ in the vicinity of $`x=0`$ are irrelevant because the approximate integral will differ but exponentially from the exact result. The integrals of the second type contain $`d\phi _1/dx`$ in the integrand. For them a step-wise approximation Eq. (23) with its derivative equal identical zero everywhere except for $`x=0`$ is an inadmissible choice. So, to keep the exponential accuracy in this case, one has to get back to Eq. (18).

The eigenvalue $`\lambda _1`$ corresponding to the approximate eigenfunction $`\phi _1`$ from Eq. (23) is evaluated via formula (15) that can be rewritten as

$$\lambda _1=_1^1\psi _0\left(\widehat{𝑱}\phi _1\right)^2=\frac{1}{R}_0^1e^{\sigma x^2}(1x^2)\left(\frac{d\phi _1}{dx}\right)^2𝑑x.$$
(24)
Substituting the derivative from Eq. (18), one finds

$$\lambda _1=C=(2/\sqrt{\pi })\sigma ^{1/2}/RG,$$

and using expression (12) for $`R`$ finally arrives at

$$\lambda _1=(4/\sqrt{\pi })\sigma ^{3/2}e^\sigma /G^2=\lambda _\text{B}/G^2.$$
(25)
With $`G`$ expanded in powers of $`\sigma ^1`$, see Eq. (12), this formula reproduces the asymptotic expression derived by Brown in Ref. Brown-P:77, . At $`G=1`$ it reduces to his initial result,Brown-PR:63  corresponding to the above-given Eq. (10). Function $`\lambda _1(\sigma )`$ from Eq. (25) is shown in Fig. 2 in comparison with the exact result obtained by a numerical solution. Indeed, at $`\sigma 3`$ the results virtually coincide.

According to expansion (5), each decrement $`\lambda _{\mathrm{}}`$ defines the reference relaxation time

$$\tau _{\mathrm{}}=2\tau _D/\lambda _{\mathrm{}}.$$
(26)
Thence from Eq. (25) we get

$$\tau _1=2\tau _D/\lambda _1=\tau _\text{B}G^2,\tau _\text{B}2\tau _D/\lambda _\text{B},$$
(27)
where $`\tau _\text{B}`$ denotes the asymptotic relaxation time obtained by Brown in Ref. Brown-PR:63, . Substituting in Eq. (27) the explicit asymptotic series (12) for $`G`$, one gets

$$\tau _1=\tau _D\frac{\sqrt{\pi }e^\sigma }{2\sigma ^{3/2}}\left(1+\frac{1}{\sigma }+\frac{7}{4\sigma ^2}+\frac{9}{2\sigma ^3}+\mathrm{}\right).$$
(28)
### II.4 Asymptotic integral time

The decrements $`\lambda _{\mathrm{}}`$ or, equivalently, relaxation times $`\tau _{\mathrm{}}`$, being the characteristics of the eigenfunctions of the distribution function, are not observable if taken as separate quantities. However, in combination they are involved in a useful directly measurable quantity, the so-called integral relaxation time. In terms of correlation functions this characteristics is defined as

$$\tau _{\mathrm{int}}=_0^{\mathrm{}}\frac{m(t)m(0)_0}{m^2(0)_0}𝑑t=_0^{\mathrm{}}\frac{x(t)x(0)_0}{x^2(0)_0}𝑑t,$$
(29)
where the angular brackets stand for the statistical ensemble averaging over the equilibrium distribution (12). As follows from Eq. (29), the integral relaxation time equals the area under the normalized decay of magnetization.

The Green function of Eq. (2), i.e., the probability density of a state $`(x,t)`$, provided the initial state is $`(x_0,0)`$, writes

$$W(x,t;x_0,0)=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\psi _{\mathrm{}}(x)\phi _{\mathrm{}}(x_0)e^{\lambda _{\mathrm{}}t}.$$
(30)
Similarly to Eq. (7), we expand the eigenfunctions in Legendre polynomials as

$$\psi _{\mathrm{}}=\frac{1}{2}\underset{k=1}{\overset{\mathrm{}}{}}(2k+1)b_k^{(\mathrm{})}P_k(x),\phi _{\mathrm{}}=\underset{k=1}{\overset{\mathrm{}}{}}a_k^{(\mathrm{})}P_k(x).$$
(31)
and introduce special notations for the first two functions

$`\psi _0=\frac{1}{2}\underset{k=0}{\overset{\mathrm{}}{}}(2k+1)S_kP_k(x),`$ (32)
$`\psi _1=\frac{1}{2}\underset{k=0}{\overset{\mathrm{}}{}}(2k+1)Q_kP_k(x).`$
The procedures to evaluate the coefficients $`S_k`$ and $`Q_k`$ and the explicit asymptotic forms for $`Q_1`$ and $`S_2`$ are given in Appendix A; note representation (11) for the equilibrium function $`\psi _0`$.

Due to Eq. (14), the coefficients in formulas (31) are related to each other by $`b_k^{(\mathrm{})}=P_kP_k^{}_0a_k^{}^{(\mathrm{})}`$. In those terms one gets for the correlator in Eq. (14):

$$x(t)x(0)_0=xx_0\psi _0W(x,t;x_0,0)𝑑x𝑑x_0=\underset{\mathrm{}=1}{}\left[b_1^{(\mathrm{})}\right]^2e^{\lambda _{\mathrm{}}t/2\tau _D},$$
(33)
where the averaging over the current coordinate $`x`$ is performed with the function $`W`$ from Eq. (30) whereas that over the initial conditions—with the equilibrium function $`\psi _0`$. Substituting expression (33) in Eq. (29) one gets the integral time in the form

$$\tau _{\mathrm{int}}=\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\tau _{\mathrm{}}\left[b_1^{(\mathrm{})}\right]^2/\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\left[b_1^{(\mathrm{})}\right]^2=\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\tau _{\mathrm{}}\left[b_1^{(\mathrm{})}\right]^2/x^2_0.$$
(34)
Unlike $`\tau _1`$, which in principle cannot be evaluatedCoff-ACP:98  analytically at arbitrary $`\sigma `$, for $`\tau _{\text{int}}`$ an exact solution is possible for arbitrary values of the anisotropy parameter. Recently two ways were proposed to obtain quadrature formulas for $`\tau _{\text{int}}`$. One methodGaIs-TMF:90  implies a direct integration of the Fokker-Planck equation. Another methodCoCr-PRE:94  involves solving three-term recurrence relations for the statistical moments of $`W`$. The emerging solution for $`\tau _{\text{int}}`$ can be expressed in a finite form in terms of hypergeometric (Kummer’s) functions. Equivalence of both approaches was proven in Ref. CoCr-PRE:96, .

In the present study, as mentioned, we are dealing in the high-barrier approximation. In this limiting case $`\lambda _1`$ is exponentially small, so that the term with $`\mathrm{}=1`$ in the numerator in Eq. (34) is far greater than the others. With allowance for Eq. (32) it can be written as

$$\tau _{\mathrm{int}}=\tau _1\left[b_1^{(1)}\right]^2/x^2_0=\tau _1Q_1^2/x^2_0.$$
(35)
The equilibrium moment calculated by definition writes

$$x^2_0=(1/2\sigma )(e^\sigma 1)=1/G1/2\sigma ,$$
(36)
and for $`\sigma 1`$, using formula (89) of Appendix A we get

$$Q_11/G.$$
(37)
Substitution of Eqs. (36) and (37) in (35) with allowance for relationships (12), (25) and (27) gives the asymptotic representation in the form

$`\tau _{\mathrm{int}}`$ $`=`$ $`\tau _\text{B}{\displaystyle \frac{2\sigma G}{(2\sigma G)}}`$
$`=`$ $`\tau _D{\displaystyle \frac{\sqrt{\pi }e^\sigma }{2\sigma ^{3/2}}}\left(1+{\displaystyle \frac{1}{\sigma }}+{\displaystyle \frac{3}{2\sigma ^2}}+{\displaystyle \frac{13}{4\sigma ^3}}+\mathrm{}\right).`$
As it is seen from formulas (28) and (II.4) written with the accuracy up to $`\sigma ^3`$, the asymptotic expressions for the interwell and integral times deviate beginning with the term $`\sigma ^2`$. This contradicts the only known to us asymptotic expansion of $`\tau _{\text{int}}`$ given in Eq. (60) of Refs. CoCr-PRE:94,  and repeated in Eq. (7.4.3.22) of the book CoKa-TLE:96 . The latter expression written with the accuracy $`\sigma ^2`$, instead of turning into Eq. (II.4) coincides with the Brown’s expression (28) for $`\tau _1`$. Meanwhile, as it follows from formula (35), such a coincidence is impossible and therefore Eqs. (60) of Refs. CoCr-PRE:94,  and (7.4.3.22) of Ref. CoKa-TLE:96,  are misleading. The necessity to rectify this issue made us to begin the demonstration of our approach with the case of the integral relaxation time. Further on we consistently apply our procedure for description of the nonlinear (third- and fifth-harmonic) dynamic susceptibilities of a solid superparamagnetic dispersion.

## III Perturbative expansions for the distribution function

### III.1 Static probing field

To find the nonlinear susceptibilities, one has to take into account the changes that the probing field induces in the basic state of the system. In the limit $`\sigma 1`$, which we deal in, the relaxation time $`\tau _1`$ of the interwell mode $`\psi _1`$ is far greater than all the other relaxation times $`\tau _k`$. This means that with respect to the intrawell modes the distribution function is in equilibrium. So it suffices to determine the effect of the probing field $`𝑯=H𝒉`$ just on $`\psi _0`$ and $`\psi _1`$. Assuming the energy function in the form

$$U+U_H=Kv(𝒆𝒏)^2IvH(𝒆𝒉),$$
(39)
\[compare with Eq. (1)\], and separating variables in Eq. (2), one arrives at the eigenfunction problem

$$\widehat{L}f_\beta =\xi \widehat{V}f_\beta ,$$
(40)
where $`\xi =IvH/kT`$ and notation $`f_\beta `$ refers to the distribution function modes that stem from $`\psi _0`$ or $`\psi _1`$ at $`H0`$, i.e., $`\beta =0`$ or 1. In Eq. (40) operator $`\widehat{L}`$ is defined by Eq. (6) whilst $`\widehat{V}=\xi \widehat{𝑱}(𝒆\times 𝒉)`$ is the operator caused by the energy term $`U_H`$ in (Eq. 39). As in above, for the non-self-conjugated spectral problem (40) we introduce the family of conjugated functions $`g_\beta `$ and set $`f_\beta =g_\beta \psi _0`$.

Following our approach, in the low-temperature limit ($`\sigma 1`$) we set to zero the eigenvalues corresponding to both $`f_0`$ and $`f_1`$; compare with Eqs. (17) and (18) for $`\psi _0`$ and $`\psi _1`$. Assuming the temperature-scaled magnetic field $`\xi `$ to be small, we treat $`U_H`$ as a perturbation Hamiltonian and expand the principal eigenfunctions as

$$f_0=\underset{n=0}{}\xi ^nf_0^{(n)},f_1=\underset{n=0}{}\xi ^nf_1^{(n)}.$$
(41)
Thence for the field-free ($`H=0`$) case one has $`f_0^{(0)}=\psi _0`$ and $`f_1^{(0)}=\psi _1`$. The same kind of expansion is assumed for $`g_\beta `$ with $`g_0^{(0)}=1`$ and $`g_1^{(0)}=\phi _1`$. Note also that in order to retain the normalizing condition we require that $`f_\beta ^{(n)}`$ have zero averages.

Substituting expansion (41) in Eq. (40) and collecting the terms of the same order in $`\xi `$, we arrive at the recurrence relation

$$\widehat{L}f_\beta ^{(n)}=\widehat{V}f_\beta ^{(n1)},$$
(42)
that for the particular cases $`\beta =0`$ and 1 with the aid of the identity $`𝒆\times 𝒉=\widehat{𝑱}(𝒆𝒉)`$ takes the formes

$`\widehat{𝑱}\psi _0\widehat{𝑱}g_0^{(n)}=\widehat{𝑱}\psi _0g_0^{(n1)}\widehat{𝑱}(𝒆𝒉),`$ (43)
$`\widehat{𝑱}\psi _0\widehat{𝑱}g_1^{(n)}=\widehat{𝑱}\psi _0g_1^{(n1)}\widehat{𝑱}(𝒆𝒉),`$
respectively. Set (43) solves easily for $`g_0`$ since $`g_0^{(0)}=\phi _0=1`$. Starting with $`n=0`$, one gets sequentially

$`g_0^{(1)}`$ $`=`$ $`(𝒆𝒉),`$ (44)
$`g_0^{(2)}`$ $`=`$ $`\frac{1}{2}\left[(𝒆𝒉)^2(𝒆𝒉)^2_0\right],`$
$`g_0^{(3)}`$ $`=`$ $`\frac{1}{6}(𝒆𝒉)^3\frac{1}{2}(𝒆𝒉)(𝒆𝒉)^2_0,`$
$`g_0^{(4)}`$ $`=`$ $`\frac{1}{24}\left[(𝒆𝒉)^4(𝒆𝒉)^4_0\right]\frac{1}{4}\left[(𝒆𝒉)^2(𝒆𝒉)^2_0(𝒆𝒉)^2_0^2\right],`$
$`g_0^{(5)}`$ $`=`$ $`\frac{1}{120}(𝒆𝒉)^5\frac{1}{12}(𝒆𝒉)^3(𝒆𝒉)^2_0\frac{1}{24}(𝒆𝒉)\left[(𝒆𝒉)^4_06(𝒆𝒉)^2_0^2\right].`$
All the obtained functions are constructed in such a way that the corresponding $`f_\beta ^{(n)}`$ satisfy the above-mentioned zero average requirement. We remark also that there is no problem to continue the calculational procedure to any order.

Evaluation of $`g_1`$ is done in two steps. At the first one, we set $`g_1^{(0)}`$ equal to the antisymmetric step-wise function (23) and its derivative equal zero. After that from the second of Eqs. (43) we can express $`g_1^{(k)}`$ in closed form. Taken up to the fourth order these “zero-derivative” solutions write

$`g_1^{(1)}`$ $`=`$ $`\phi _1(𝒆𝒉)\phi _1(𝒆𝒉)_0,`$ (45)
$`g_1^{(2)}`$ $`=`$ $`\frac{1}{2}\phi _1(𝒆𝒉)^2(𝒆𝒉)\phi _1(𝒆𝒉)_0,`$
$`g_1^{(3)}`$ $`=`$ $`\frac{1}{6}\left[\phi _1(𝒆𝒉)^3\phi _1(𝒆𝒉)^3_0\right]\frac{1}{2}\phi _1(𝒆𝒉)_0\left[(𝒆𝒉)^2(𝒆𝒉)^2_0\right],`$
$`g_1^{(4)}`$ $`=`$ $`\frac{1}{24}\phi _1(𝒆𝒉)^4\frac{1}{6}\phi _1(𝒆𝒉)_0\left[(𝒆𝒉)^33(𝒆𝒉)(𝒆𝒉)^2_0\right]`$
$`\frac{1}{6}(𝒆𝒉)\phi _1(𝒆𝒉)^3_0.`$
Note the alternating parity in $`𝒆`$ with the term order growth in both Eqs. (44) and (45).

It is instructive to compare the approximate expressions (45) with the numerical results obtained without simplification of $`g_1^{(0)}`$. To be specific, we consider the case when probing field is applied along the particle easy axis $`𝒏`$. Then Eqs. (43) become one-dimensional and the second of them writes

$$\frac{dg_1^{(n)}}{dx}=g_1^{(n1)}.$$
(46)
Its “zero-derivative” solutions up to the second order follow from the first two lines of Eqs. (45):

$$g_1^{(1)}=\phi _1x\phi _1x_0,g_1^{(2)}=\frac{1}{2}\phi _1x^2x\phi _1x_0.$$
(47)
In Figs. 3 and 4 these functions are compared to the numerical solutions of Eq. (46). For our calculation, the most important is the behavior of those functions near $`x=\pm 1`$ since these regions yield the main contribution when integrated with the weight function $`\psi _0`$. As one can see from the figures, the “zero-derivative” solution $`g_1^{(1)}`$ agrees well with the exact one, whilst $`g_1^{(2)}`$ deviates significantly. This discrepancy is due to the change of the barrier height that occurs in the second order with respect to the probing field amplitude, and manifests itself in all the even orders of the perturbation expansion. Correction of solution (47) makes the second step of our procedure. For that we integrate Eq. (46) two times by parts

and substitute there the “zero-derivative” form of $`g_1^{(1)}`$ from Eq. (47):

$$g_1^{(2)}=\frac{1}{2}x^2\phi _1xx\phi _1+\frac{1}{2}x^2\frac{d\phi _1}{dx}𝑑x,$$
(48)
Thus one finds that the corrected $`g_1^{(2)}`$ differs from this of Eq. (47) by adding a step-wise \[alike that of Eq. (23)\] term

$$g_1^{(2)}=\frac{1}{2}x^2\phi _1xx\phi _1+D_2\phi _1,$$
(49)
with the amplitude

$$D_2=\frac{1}{2}_0^1x^2\frac{d\phi _1}{dx}𝑑x.$$
(50)
We remark that the results of evaluation of the integrals $`I_{2k}=_0^1x^{2k}(d\phi _1/dx)𝑑x`$ can be arranged in the table

$$\begin{array}{cccc}k& \mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}& 1& 2\\ & & & \\ I_{2k}& \mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}& \mathrm{\hspace{0.25em}1}G^1& \mathrm{\hspace{0.25em}1}(1+1/2\sigma )G^1\end{array}.$$
(51)
so that Eq. (50) gives

$$D_2=\frac{1}{2}I_{2k}=\frac{G1}{2G}=\frac{1}{4\sigma }+\frac{1}{4\sigma ^2}+\frac{5}{8\sigma ^3}+\frac{37}{16\sigma ^4}+\mathrm{}$$
(52)
Function $`g_1^{(2)}`$ corrected in such a way is shown in Fig. 4 by asterisks. It is seen that the corrected dependence with a fairly good accuracy follows the the numerically obtained curve.

In a similar way one can prove that the corrected function $`g_1^{(4)}`$ has the form

$`g_1^{(4)}={\displaystyle \frac{1}{24}}\phi _1x^4{\displaystyle \frac{1}{6}}\left[\phi _1x_0x^33xx^2_0\right]`$ (53)
$`{\displaystyle \frac{1}{6}}x\phi _1x^3_0+D_2g_1^{(2)}+D_4\phi _1,`$
where the corrected function $`g_1^{(2)}`$ given by Eq. (49) is used and

$`D_4={\displaystyle \frac{1}{24}}I_4D_2^2={\displaystyle \frac{10\sigma G^222\sigma G+G+12\sigma }{48\sigma G^2}}`$ (54)
$`={\displaystyle \frac{1}{32\sigma ^2}}{\displaystyle \frac{1}{16\sigma ^3}}{\displaystyle \frac{5}{32\sigma ^4}}{\displaystyle \frac{29}{64\sigma ^5}}+\mathrm{}`$
In the general case, when the direction of the probing field does not coincide with the particle anisotropy axis, the corrected functions $`g_1^{(n)}`$ still can be written as

$`g_1^{(2)}`$ $`=`$ $`\frac{1}{2}\phi _1(𝒆𝒉)^2(𝒆𝒉)\phi _1(𝒆𝒉)_0+D_2\phi _1,`$
$`g_1^{(3)}`$ $`=`$ $`\frac{1}{6}\left[\phi _1(𝒆𝒉)^3\phi _1(𝒆𝒉)^3_0\right]\frac{1}{2}\phi _1(𝒆𝒉)_0\left[(𝒆𝒉)^2(𝒆𝒉)^2_0\right]+D_2g_1^{(1)},`$
$`g_1^{(4)}`$ $`=`$ $`\frac{1}{24}\phi _1(𝒆𝒉)^4\frac{1}{6}\phi _1(𝒆𝒉)_0\left[(𝒆𝒉)^33(𝒆𝒉)(𝒆𝒉)^2_0\right]`$ (55)
$`\frac{1}{6}(𝒆𝒉)\phi _1(𝒆𝒉)^3_0+D_2g_1^{(2)}+D_4\phi _1.`$
But since Eqs. (43) cannot be reduced to a form like Eq. (46), the correcting coefficients $`D_2`$ and $`D_4`$ cannot be presented in a closed form. In this case the corrected solutions taking into account the behavior of function $`\phi _1`$ around zero are built up as power series near $`x=0`$; such a procedure for the coefficients $`D_2`$ and $`D_4`$ is described in Appendix B.

### III.2 Dynamic probing field

To obtain the dynamic susceptibilities, one has to find the distribution function $`W`$ in the oscillating probing field $`\xi \mathrm{exp}(i\omega t)`$. For this situation the kinetic equation (2) takes the form

$$\left(2\tau _D\frac{}{t}+\widehat{L}\right)W(t)=\xi \widehat{V}e^{i\omega t}W(t),$$
(56)
where the operators $`\widehat{L}`$ and $`\widehat{V}`$ have been introduced in above. Assuming that the exciting field amplitude is not too high, we expand the steady-state oscillatory solution of Eq. (56) in a power series with respect to $`\xi `$:

$$W(t)=\underset{n=0}{}\xi ^nW^{(n)}e^{in\omega t}.$$
(57)
Note that, mathematically, representation (57) is not complete. Indeed, in a general case the exact amplitude of the $`n\omega `$-mode must contain, along with the contribution $`\xi ^n`$, an infinite set of terms $`\xi ^{n+2}`$, $`\xi ^{n+4}`$, etc. However, in a weak field limit $`\xi <1`$ the terms with higher powers are of minor importance so that the main contribution to the magnetization response signal filtered at the frequency $`n\omega `$ is proportional to $`\xi ^n`$.

Substituting Eq. (57) in (56) we arrive at the recurrence set

$$\left(2in\omega \tau _D+\widehat{L}\right)W^{(n)}=\widehat{V}W^{(n1)},$$
(58)
that we solve sequentially starting from $`n=1`$. At the first step the function in the right-hand side corresponds to the equilibrium case ($`\xi =0`$). Therefore, $`W^{(0)}=\psi _0`$, where the latter function is defined by Eq. (11) and is frequency-independent. Combining Eq. (42) written down for $`\beta =0`$ and $`n=1`$ and Eq. (58), we eliminate the operator $`\widehat{V}`$ and get

$$\left(2i\omega \tau _D+\widehat{L}\right)W^{(1)}=\widehat{L}f_0^{(1)}.$$
(59)
Now we expand the functions subjected to operator $`\widehat{L}`$ with respect to the set $`\{\psi _k\}`$ of its eigenfunctions, see Eq. (6):

$$W^{(1)}=c_j^{(1)}(\omega )\psi _jf_0^{(1)}=(\phi _j\left|f_0^{(1)}\right)\psi _j;$$
(60)
here $`(\phi |f)`$ denotes functional scalar multiplication, i.e., the integral of the product $`\phi f`$ over all the orientations of $`𝒆`$. Substitution of Eq. (60) in (59), multiplication of it from the left by $`\phi _k`$ and integration, render the expansion coefficient as

$$c_k^{(1)}(\omega )=(\phi _k\left|f_0^{(1)}\right)[1+i\omega \tau _k]^1,$$
(61)
where the reference relaxation times are defined by Eq. (26).

In the low-frequency limit only $`\omega \tau _1`$ is set to be non-zero whilst all the higher modes are taken at equilibrium ($`\omega \tau _k=0`$). Thence, when constructing $`W^{(1)}`$ via Eq. (60), by adding and subtracting a term with $`c_1^{(1)}(0)`$, one can present the first-order solution in the form

$$W^{(1)}=f_0^{(1)}\frac{i\omega \tau _1}{1+i\omega \tau _1}(\phi _1\left|f_0^{(1)}\right)\psi _1,$$
(62)
where $`f_0^{(1)}`$, as seen from Eq. (59), is the equilibrium solution for the same value of the field amplitude $`\xi `$. We remind that the functions without upper index belong to the fundamental set defined by Eqs. (6) whereas those with an upper index are evaluated in the framework of the perturbation scheme described in Sec. III.1.

In the next order in $`\xi `$ the function $`W^{(1)}`$ is substituted in the right-hand side of Eq. (58) and through a procedure alike to that leading to Eqs. (59)–(61), the function $`W^{(2)}`$ is found. We carry on this cycle up to $`k=5`$. The results write

$$W^{(2)}=f_0^{(2)}\frac{i\omega \tau _1}{1+i\omega \tau _1}(\phi _1\left|f_0^{(1)}\right)f_1^{(1)},$$
(63)
$`W^{(3)}=f_0^{(3)}+[(\phi _1\left|f_0^{(1)}\right)(\phi _1\left|f_1^{(2)}\right)\left(\phi _1\left|f_0^{(3)}\right)\right]f_1^{(0)}(\phi _1\left|f_0^{(1)}\right)f_1^{(2)}`$
$`+{\displaystyle \frac{1}{1+i\omega \tau _1}}[(\phi _1\left|f_0^{(1)}\right)f_1^{(2)}\frac{3}{2}(\phi _1\left|f_0^{(1)}\right)\left(\phi _1\left|f_1^{(2)}\right)f_1^{(0)}\right]`$
$`+{\displaystyle \frac{1}{1+3i\omega \tau _1}}[(\phi _1\left|f_0^{(3)}\right)+\frac{1}{2}(\phi _1\left|f_0^{(1)}\right)\left(\phi _1\left|f_1^{(2)}\right)\right]f_1^{(0)},`$ (64)
$`W^{(4)}=f_0^{(4)}+[(\phi _1\left|f_0^{(1)}\right)(\phi _1\left|f_1^{(2)}\right)\left(\phi _1\left|f_0^{(3)}\right)\right]f_1^{(1)}(\phi _1\left|f_0^{(1)}\right)f_1^{(3)}`$
$`+{\displaystyle \frac{1}{1+i\omega \tau _1}}[(\phi _1\left|f_0^{(1)}\right)f_1^{(3)}\frac{3}{2}(\phi _1\left|f_0^{(1)}\right)\left(\phi _1\left|f_1^{(2)}\right)f_1^{(1)}\right]`$
$`+{\displaystyle \frac{1}{1+3i\omega \tau _1}}[(\phi _1\left|f_0^{(3)}\right)+\frac{1}{2}(\phi _1\left|f_0^{(1)}\right)\left(\phi _1\left|f_1^{(2)}\right)\right]f_1^{(1)},`$ (65)
$`W^{(5)}`$ $`=`$ $`f_0^{(5)}(\phi _1\left|f_0^{(5)}\right)f_1^{(0)}+(\phi _1\left|f_0^{(1)}\right)[(\phi _1\left|f_1^{(4)}\right)f_1^{(0)}f_1^{(4)}]`$ (66)
$`+`$ $`[(\phi _1\left|f_0^{(1)}\right)(\phi _1\left|f_1^{(2)}\right)\left(\phi _1\left|f_0^{(3)}\right)\right][f_1^{(2)}\left(\phi _1\left|f_1^{(2)}\right)f_1^{(0)}\right]`$
$`+`$ $`{\displaystyle \frac{1}{1+i\omega \tau _1}}(\phi _1\left|f_0^{(1)}\right)\{f_1^{(4)}\frac{3}{2}(\phi _1\left|f_1^{(2)}\right)f_1^{(2)}+[\frac{15}{8}\left(\phi _1\left|f_1^{(2)}\right)^2\frac{5}{4}\left(\phi _1\left|f_1^{(4)}\right)\right]f_1^{(0)}\right\}`$
$`+`$ $`{\displaystyle \frac{1}{1+3i\omega \tau _1}}[(\phi _1\left|f_0^{(3)}\right)+\frac{1}{2}(\phi _1\left|f_0^{(1)}\right)\left(\phi _1\left|f_1^{(2)}\right)\right][f_1^{(2)}\frac{5}{2}\left(\phi _1\left|f_1^{(2)}\right)f_1^{(0)}\right]`$
$`+`$ $`{\displaystyle \frac{1}{1+5i\omega \tau _1}}[(\phi _1\left|f_0^{(5)}\right)+\frac{1}{4}(\phi _1\left|f_0^{(1)}\right)(\phi _1\left|f_1^{(4)}\right)+\frac{3}{8}(\phi _1\left|f_0^{(1)}\right)(\phi _1\left|f_1^{(2)}\right)^2`$
$`+\frac{3}{2}(\phi _1\left|f_0^{(3)}\right)\left(\phi _1\left|f_1^{(2)}\right)\right]f_1^{(0)}.`$
We remark an important feature of Eqs. (63)–(66): they do not contain dispersion factors of even orders. This ensures that the frequency dependence of the full distribution function $`W`$ incorporates only dispersion factors with odd multiples of the basic frequency. Qualitatively, this is the result of absence of the interwell mode for the statistical moments of even orders. Technically, it is due to vanishing of the products $`(\phi _1\left|f_k^{(\mathrm{})}\right)`$ entering Eqs. (62)–(66) if the sum $`k+\mathrm{}`$ is even. This rule follows immediately from combination of the oddity of $`\phi _1`$, see Sec. II, with the parity properties of the functions $`f_k^{(\mathrm{})}`$ introduced in Sec. III.1.

For actual calculations one needs the values of the scalar products entering Eqs. (62)–(66). In Appendix C we obtain their representations in terms of the moments $`Q_k`$ and $`S_k`$ of the functions $`\psi _0`$ and $`\psi _1`$, respectively. The procedures of asymptotic expansion of $`Q_k`$ and $`S_k`$ are given in Appendix A.

## IV Dynamic susceptibilities

The set of magnetic susceptibilities of an assembly of non-interacting particles with the number density $`c`$ is defined by the relation

$$M=\chi ^{(1)}H+\chi ^{(3)}H^3+\chi ^{(5)}H^5+\mathrm{}$$
(67)
that describes the magnetization of the system in the direction of the probing field $`𝑯=H𝒉`$. Therefore, of all the components of the corresponding susceptibility tensors, we retain the combinations that determine the response in the direction of the probing field. With representation (57) for the distribution function, this magnetization component takes the form

$`M=cIv(𝒆𝒉)`$ (68)
$`=c{\displaystyle \underset{n=1}{}}H^n{\displaystyle \frac{I^{n+1}v^{n+1}}{(kT)^n}}e^{in\omega t}{\displaystyle (𝒆𝒉)W^{(n)}𝑑𝒆},`$
and the susceptibilities can be found by a direct comparison with Eq. (67). In other words, the set of $`\chi ^{(n)}`$ is expressed through the perturbation functions $`W^{(n)}`$ found in the preceding section. Therefore, evaluation of $`\chi ^{(n)}`$ becomes, although tedious, but simple procedure. Remarkably, the final expressions come out rather compact.

### IV.1 Linear susceptibility

The resulting expression can be presented in the form

$$\chi _\omega ^{(1)}=\chi _0^{(1)}\left(B_0^{(1)}+\frac{B_1^{(1)}}{1+i\omega \tau _1}\right),\chi _0^{(1)}=\frac{cI^2v^2}{3kT},$$
(69)
which follows from substituting Eq. (62) in (68). Each of the two frequency-independent coefficients $`B^{(1)}`$, being the result of statistical averaging over the orientational variable $`𝒆`$, see Appendix C, expands into a series of Legendre polynomials with respect to $`\beta `$, the angle between the direction $`𝒉`$ of the probing field and the particle easy axis $`𝒏`$. This can be written as

$`\{\begin{array}{cc}B_0^{(1)}=b_{00}^{(1)}+b_{02}^{(1)}P_2(\mathrm{cos}\beta ),\hfill & \\ B_1^{(1)}=b_{10}^{(1)}+b_{12}^{(1)}P_2(\mathrm{cos}\beta ),\hfill & \end{array}`$ (70)
$`\left(\begin{array}{cc}b_{00}^{(1)}& b_{02}^{(1)}\\ b_{10}^{(1)}& b_{12}^{(1)}\end{array}\right)=\left(\begin{array}{cc}1Q_1^2& 2S_22Q_1^2\\ Q_1^2& 2Q_1^2\end{array}\right).`$
Definitions of functions $`Q_1`$ and $`S_2`$ and their explicit asymptotic representations are given in Appendix A. The asymptotic series for the coefficients $`b_{\alpha \beta }^{(1)}`$ derived on the base of expansion (12) and Eq. (37) are

$`b_{02}^{(1)}={\displaystyle \frac{1}{\sigma }}+{\displaystyle \frac{1}{4\sigma ^3}}+{\displaystyle \frac{13}{8\sigma ^4}}+{\displaystyle \frac{165}{16\sigma ^5}}+{\displaystyle \frac{2273}{32\sigma ^6}}+{\displaystyle \frac{34577}{64\sigma ^7}}+{\displaystyle \frac{581133}{128\sigma ^8}}+\mathrm{}`$ (75)
$`b_{10}^{(1)}=1{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{3}{4\sigma ^2}}{\displaystyle \frac{2}{\sigma ^3}}{\displaystyle \frac{31}{4\sigma ^4}}{\displaystyle \frac{153}{4\sigma ^5}}{\displaystyle \frac{3629}{16\sigma ^6}}{\displaystyle \frac{1564}{\sigma ^7}}{\displaystyle \frac{785931}{64\sigma ^8}}+\mathrm{}`$
The other components, namely, $`b_{00}^{(1)}`$ and $`b_{12}^{(1)}`$, may be constructed straightforwardly using their relations with the given ones, see Eqs. (70). For a random system, that is for an assembly of non-interacting particles with a chaotic distribution of the anisotropy axes, the average of any Legendre polynomial is zero, so that $`\overline{B^{(1)}}_k=b_{k0}^{(1)}`$, and the linear dynamic susceptibility reduces to

$$\overline{\chi ^{(1)}}_\omega =\chi _0^{(1)}\frac{1+i\omega \tau _1b_{00}}{1+i\omega \tau _1},$$
(76)
that is the asymptotic representation of the full expression given by formula (39) of Ref. RaSt-PRB:97 .

### IV.2 Cubic susceptibility

As follows from definitions (67) and (68), the third-order susceptibility is defined through the response at the triple frequency that at weak $`H`$ scales as $`H^3`$. Performing calculations along the same scheme as for $`\chi ^{(1)}`$, one arrives at the sum-of-relaxators representation

$`\chi _{3\omega }^{(3)}={\displaystyle \frac{1}{4}}\chi _0^{(3)}\left(B_0^{(3)}+{\displaystyle \frac{B_1^{(3)}}{1+i\omega \tau _1}}+{\displaystyle \frac{B_3^{(3)}}{1+3i\omega \tau _1}}\right),`$ (77)
$`\chi _0^{(3)}={\displaystyle \frac{cI^4v^4}{(kT)^3}},`$
where the coefficients expand as

$$B_k^{(3)}=b_{k0}^{(3)}+b_{k2}^{(3)}P_2(\mathrm{cos}\beta )+b_{k4}^{(3)}P_4(\mathrm{cos}\beta ),k=0,1,3\mathrm{}$$
(78)
up to the fourth Legendre polynomial in $`\mathrm{cos}\beta `$.

The explicit expansions for the amplitudes $`b_{\alpha \beta }^{(3)}`$ are

$`b_{00}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{1}{30\sigma ^3}}+{\displaystyle \frac{47}{240\sigma ^4}}+{\displaystyle \frac{49}{40\sigma ^5}}+{\displaystyle \frac{815}{96\sigma ^6}}+{\displaystyle \frac{7837}{120\sigma ^7}}+{\displaystyle \frac{355391}{640\sigma ^8}}+\mathrm{}`$ (79)
$`b_{02}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{1}{42\sigma ^3}}+{\displaystyle \frac{2}{21\sigma ^4}}+{\displaystyle \frac{4}{7\sigma ^5}}+{\displaystyle \frac{1385}{336\sigma ^6}}+{\displaystyle \frac{11231}{336\sigma ^7}}+{\displaystyle \frac{19083}{64\sigma ^8}}+\mathrm{}`$
$`b_{04}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{2}{35\sigma ^3}}{\displaystyle \frac{8}{35\sigma ^4}}{\displaystyle \frac{41}{35\sigma ^5}}{\displaystyle \frac{50}{7\sigma ^6}}{\displaystyle \frac{1756}{35\sigma ^7}}{\displaystyle \frac{63749}{160\sigma ^8}}+\mathrm{}`$
$`b_{10}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{1}{15}}{\displaystyle \frac{1}{6\sigma }}{\displaystyle \frac{23}{240\sigma ^2}}{\displaystyle \frac{61}{192\sigma ^3}}{\displaystyle \frac{1357}{960\sigma ^4}}{\displaystyle \frac{235447}{30720\sigma ^5}}{\displaystyle \frac{11962691}{245760\sigma ^6}}{\displaystyle \frac{694849241}{1966080\sigma ^7}}{\displaystyle \frac{15133953221}{5242880\sigma ^8}}+\mathrm{}`$
$`b_{12}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{13}{84}}{\displaystyle \frac{65}{168\sigma }}{\displaystyle \frac{25}{168\sigma ^2}}{\displaystyle \frac{863}{1344\sigma ^3}}{\displaystyle \frac{3931}{1344\sigma ^4}}{\displaystyle \frac{698911}{43008\sigma ^5}}{\displaystyle \frac{35309123}{344064\sigma ^6}}{\displaystyle \frac{2061480665}{2752512\sigma ^7}}{\displaystyle \frac{45071465669}{7340032\sigma ^8}}+\mathrm{}`$
$`b_{14}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{1}{35}}{\displaystyle \frac{1}{14\sigma }}+{\displaystyle \frac{2}{35\sigma ^2}}{\displaystyle \frac{1}{112\sigma ^3}}{\displaystyle \frac{73}{560\sigma ^4}}{\displaystyle \frac{17033}{17920\sigma ^5}}{\displaystyle \frac{1007549}{143360\sigma ^6}}{\displaystyle \frac{64390439}{1146880\sigma ^7}}{\displaystyle \frac{4493994417}{9175040\sigma ^8}}+\mathrm{}`$
$`b_{30}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{2}{15}}+{\displaystyle \frac{3}{10\sigma }}+{\displaystyle \frac{1}{16\sigma ^2}}+{\displaystyle \frac{337}{960\sigma ^3}}+{\displaystyle \frac{499}{320\sigma ^4}}+{\displaystyle \frac{85309}{10240\sigma ^5}}+{\displaystyle \frac{2563751}{49152\sigma ^6}}+{\displaystyle \frac{245269747}{655360\sigma ^7}}+{\displaystyle \frac{47628510799}{15728640\sigma ^8}}+\mathrm{}`$
$`b_{32}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{29}{84}}+{\displaystyle \frac{43}{56\sigma }}+{\displaystyle \frac{11}{56\sigma ^2}}+{\displaystyle \frac{1279}{1344\sigma ^3}}+{\displaystyle \frac{1881}{448\sigma ^4}}+{\displaystyle \frac{320765}{14336\sigma ^5}}+{\displaystyle \frac{48133699}{34406\sigma ^6}}+{\displaystyle \frac{920146163}{91750\sigma ^7}}+{\displaystyle \frac{178560431695}{22020096\sigma ^8}}+\mathrm{}`$
$`b_{34}^{\left(3\right)}`$ $`=`$ $`{\displaystyle \frac{11}{105}}+{\displaystyle \frac{47}{210\sigma }}+{\displaystyle \frac{2}{21\sigma ^2}}+{\displaystyle \frac{559}{1680\sigma ^3}}+{\displaystyle \frac{2419}{1680\sigma ^4}}+{\displaystyle \frac{409499}{53760\sigma ^5}}+{\displaystyle \frac{4080395}{86016\sigma ^6}}+{\displaystyle \frac{1166954357}{3440640\sigma ^7}}+{\displaystyle \frac{75334335763}{27525120\sigma ^8}}+\mathrm{}`$
For a random system, the averages of Legendre polynomials drop out and $`\overline{B_k^{(3)}}=b_{k0}^{(3)}`$. With respect to formalism constructed in Ref. RaSt-PRB:97 , the above expressions yield the asymptotic representations for formulas (42) and (43) there.

### IV.3 Fifth-order susceptibility

The susceptibility of the fifth order writes in an expectable way as a sum of three relaxators:

$$\chi _{5\omega }^{(5)}=\frac{1}{16}\chi _0^{(5)}\left(B_0^{(5)}+\frac{B_1^{(5)}}{1+i\omega \tau }+\frac{B_3^{(5)}}{1+3i\omega \tau }+\frac{B_5^{(5)}}{1+5i\omega \tau }\right),\chi _0^{(5)}=\frac{cI^6v^6}{(kT)^5},$$
(80)
with the coefficients

$$B_k^{(5)}=b_{k0}^{(5)}+b_{k2}^{(5)}P_2(\mathrm{cos}\beta )+b_{k4}^{(5)}P_4(\mathrm{cos}\beta )+b_{k6}^{(5)}P_6(\mathrm{cos}\beta )k=0,1,3,5.$$
(81)
The explicit asymptotic series are

$`b_{00}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1}{80\sigma ^5}}+{\displaystyle \frac{367}{2240\sigma ^6}}+{\displaystyle \frac{123}{70\sigma ^7}}+{\displaystyle \frac{41233}{2240\sigma ^8}}+\mathrm{}`$
$`b_{10}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1}{96}}{\displaystyle \frac{19}{420\sigma }}+{\displaystyle \frac{1}{120\sigma ^2}}{\displaystyle \frac{65}{1792\sigma ^3}}{\displaystyle \frac{79}{336\sigma ^4}}{\displaystyle \frac{85913}{57344\sigma ^5}}{\displaystyle \frac{72636131}{6881280\sigma ^6}}{\displaystyle \frac{4543038053}{55050240\sigma ^7}}{\displaystyle \frac{14938598691}{20971520\sigma ^8}}+\mathrm{}`$
$`b_{30}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{47}{560}}+{\displaystyle \frac{11}{35\sigma }}{\displaystyle \frac{29}{280\sigma ^2}}+{\displaystyle \frac{437}{1920\sigma ^3}}+{\displaystyle \frac{5473}{4480\sigma ^4}}+{\displaystyle \frac{1046209}{143360\sigma ^5}}+{\displaystyle \frac{169435283}{3440640\sigma ^6}}+{\displaystyle \frac{684614895}{1835008\sigma ^7}}+{\displaystyle \frac{230861266333}{73400320\sigma ^8}}+\mathrm{}`$
$`b_{50}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{311}{3360}}{\displaystyle \frac{137}{420\sigma }}+{\displaystyle \frac{13}{105\sigma ^2}}{\displaystyle \frac{5911}{26880\sigma ^3}}{\displaystyle \frac{2141}{1920\sigma ^4}}{\displaystyle \frac{1874309}{286720\sigma ^5}}{\displaystyle \frac{299470403}{6881280\sigma ^6}}{\displaystyle \frac{17964831133}{55050240\sigma ^7}}{\displaystyle \frac{400677748549}{146800640\sigma ^8}}+\mathrm{}`$
$`b_{02}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1}{112\sigma ^5}}+{\displaystyle \frac{3}{28\sigma ^6}}+{\displaystyle \frac{507}{448\sigma ^7}}+{\displaystyle \frac{5377}{448\sigma ^8}}+\mathrm{}`$
$`b_{12}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{13}{504}}{\displaystyle \frac{19}{168\sigma }}+{\displaystyle \frac{23}{672\sigma ^2}}{\displaystyle \frac{737}{8064\sigma ^3}}{\displaystyle \frac{2959}{5376\sigma ^4}}{\displaystyle \frac{99733}{28672\sigma ^5}}{\displaystyle \frac{50499149}{2064384\sigma ^6}}{\displaystyle \frac{350973527}{1835008\sigma ^7}}{\displaystyle \frac{72765921299}{44040192\sigma ^8}}+\mathrm{}`$
$`b_{32}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{5}{21}}+{\displaystyle \frac{149}{168\sigma }}{\displaystyle \frac{193}{672\sigma ^2}}+{\displaystyle \frac{5245}{8064\sigma ^3}}+{\displaystyle \frac{18677}{5376\sigma ^4}}+{\displaystyle \frac{1785635}{86016\sigma ^5}}+{\displaystyle \frac{289305193}{2064384\sigma ^6}}+{\displaystyle \frac{5846947361}{5505024\sigma ^7}}+{\displaystyle \frac{394448762615}{44040192\sigma ^8}}+\mathrm{}`$
$`b_{52}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{139}{504}}{\displaystyle \frac{27}{28\sigma }}+{\displaystyle \frac{109}{336\sigma ^2}}{\displaystyle \frac{1343}{2016\sigma ^3}}{\displaystyle \frac{9203}{2688\sigma ^4}}{\displaystyle \frac{431321}{21504\sigma ^5}}{\displaystyle \frac{9839105}{73728\sigma ^6}}{\displaystyle \frac{196654913}{196608\sigma ^7}}{\displaystyle \frac{30690812563}{3670016\sigma ^8}}+\mathrm{}`$
$`b_{04}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{3}{140\sigma ^5}}{\displaystyle \frac{1563}{6160\sigma ^6}}{\displaystyle \frac{7767}{3080\sigma ^7}}{\displaystyle \frac{613353}{24640\sigma ^8}}+\mathrm{}`$
$`b_{14}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{15}{2464}}{\displaystyle \frac{183}{6160\sigma }}+{\displaystyle \frac{713}{24640\sigma ^2}}{\displaystyle \frac{433}{19712\sigma ^3}}{\displaystyle \frac{409}{4928\sigma ^4}}{\displaystyle \frac{319665}{630784\sigma ^5}}{\displaystyle \frac{8222083}{2293760\sigma ^6}}{\displaystyle \frac{5744848239}{201850880\sigma ^7}}{\displaystyle \frac{403943151013}{1614807040\sigma ^8}}+\mathrm{}`$
$`b_{34}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{29}{280}}+{\displaystyle \frac{293}{770\sigma }}{\displaystyle \frac{47}{385\sigma ^2}}+{\displaystyle \frac{7081}{24640\sigma ^3}}+{\displaystyle \frac{74647}{49280\sigma ^4}}+{\displaystyle \frac{7137293}{788480\sigma ^5}}+{\displaystyle \frac{385804437}{6307840\sigma ^6}}+{\displaystyle \frac{4682760003}{10092544\sigma ^7}}+{\displaystyle \frac{1580817298041}{403701760\sigma ^8}}+\mathrm{}`$
$`b_{54}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1713}{12320}}{\displaystyle \frac{2929}{6160\sigma }}+{\displaystyle \frac{2551}{24640\sigma ^2}}{\displaystyle \frac{34863}{98560\sigma ^3}}{\displaystyle \frac{92061}{49280\sigma ^4}}{\displaystyle \frac{34432191}{3153920\sigma ^5}}{\displaystyle \frac{23756287}{327680\sigma ^6}}{\displaystyle \frac{15647080587}{28835840\sigma ^7}}{\displaystyle \frac{7317549380671}{1614807040\sigma ^8}}+\mathrm{}`$
$`b_{06}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1}{616\sigma ^6}}{\displaystyle \frac{3}{77\sigma ^7}}{\displaystyle \frac{1467}{2464\sigma ^8}}+\mathrm{}`$
$`b_{16}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1}{1584}}+{\displaystyle \frac{1}{1848\sigma }}+{\displaystyle \frac{7}{1056\sigma ^2}}+{\displaystyle \frac{337}{88704\sigma ^3}}+{\displaystyle \frac{53}{1848\sigma ^4}}+{\displaystyle \frac{51433}{315392\sigma ^5}}+{\displaystyle \frac{2188103}{2064384\sigma ^6}}+{\displaystyle \frac{471762913}{60555264\sigma ^7}}+{\displaystyle \frac{4428495037}{69206016\sigma ^8}}+\mathrm{}`$
$`b_{36}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1}{84}}+{\displaystyle \frac{10}{231\sigma }}{\displaystyle \frac{17}{924\sigma ^2}}+{\displaystyle \frac{103}{3168\sigma ^3}}+{\displaystyle \frac{2489}{14784\sigma ^4}}+{\displaystyle \frac{236615}{236544\sigma ^5}}+{\displaystyle \frac{38344237}{5677056\sigma ^6}}+{\displaystyle \frac{776232845}{15138816\sigma ^7}}+{\displaystyle \frac{52467158027}{121110528\sigma ^8}}+\mathrm{}`$
$`b_{56}^{\left(5\right)}`$ $`=`$ $`{\displaystyle \frac{1207}{55440}}{\displaystyle \frac{661}{9240\sigma }}+{\displaystyle \frac{17}{12320\sigma ^2}}{\displaystyle \frac{5525}{88704\sigma ^3}}{\displaystyle \frac{1169}{3520\sigma ^4}}{\displaystyle \frac{9116467}{4730880\sigma ^5}}{\displaystyle \frac{131486063}{10321920\sigma ^6}}{\displaystyle \frac{1918435847}{20185088\sigma ^7}}{\displaystyle \frac{639291980689}{807403520\sigma ^8}}+\mathrm{},`$
(82)
and for a random system, as for the lower orders, $`\overline{B_k^{(5)}}=b_{k0}^{(0)}`$.

## V Discussion

The above derived formulas despite their hefty look are very practical. Indeed, they present the nonlinear initial susceptibilities of a superparamagnetic particulate medium as analytical expressions of arbitrary accuracy. With respect to the frequency dependence they give the exact full structure of the susceptibility and prove that it is very simple thus putting former intuitive considerations on a solid ground. This makes our formulas a handy tool for asymptotic analysis. Yet more convenient they are for numerical work because with their use the difficult and time-consuming procedure of solving the differential equations is replaced by a plain summation of certain power series. For example, if to employ Eqs. (76)–(IV.3), a computer code that fits simultaneously experimental data on linear and a set of nonlinear susceptibilities taking into account the particle polydispersity of any kind (easy axes directions, activation volume, anisotropy constants) becomes a very fast procedure.

Graphic examples justifying our claims are presented in Figs. 5 and 6, where the components of two nonlinear complex susceptibilities are plotted as the functions of the parameter $`\sigma `$. For a given sample, $`\sigma `$ in a natural way serves as a dimensionless inverse temperature. In those figures, the solid lines correspond to the above-proposed asymptotic formulas where we keep the terms up to $`\sigma ^3`$. The circles show the results of numerically-exact solutions obtained by the method described in Ref. RaSt-PRB:97, . Note that even at $`\sigma 5`$ the accuracy is still rather high.

The model that may be called the predecessor of the afore-derived results was proposed in Ref. KlYa-JMMM:98, . There, the authors calculated the initial susceptibilities up to the seventh order having replaced a superparamagnetic assembly by a two-level macrospin system. The interrelation between the present work and Ref. KlYa-JMMM:98,  closely resembles the situation with the evaluation of the rate of a superparamagnetic process. First in 1949 NéelNeel-AGCR:49  and then, ten years later, BrownBrown-JAP:59  had evaluated the superparamagnetic time in the framework of a two-level model. In such a framework, one allows for the magnetic moment flips but totally neglects its possible diffusion over energetically less-favorable directions. In 1963 BrownBrown-PR:63  had succeeded to overcome this artificial assumption and took into account the possibility for the magnetic moment to wander over all $`4\pi `$ radians.

In the present case, the obtained $`v/T`$ dependencies of the nonlinear susceptibilities and those from Ref. KlYa-JMMM:98,  are qualitatively the same. Their most typical feature is the double-peak shape. Quantitatively, however, the corresponding lines differ and do not reduce to one another in any case. Indeed, as long as the temperature is finite (whatever low), the configurational space for the unit vector $`𝒆`$ of the magnetic moment is the full ($`4\pi `$-radian) solid angle; its reduction to just two directions along a bidirectional axis could not be done otherwise than “by hand”. This is exactly what the two-level Ising-like model does: it forcibly imparts a quantum property (discrete spin projections) to a macrospin assembly. From the calculational viewpoint, another essential demerit of the resultsKlYa-JMMM:98  is that the coefficients in the susceptibility formulas are not given in an analytic form. The authors propose to evaluate them by solving an infinite set of recurrence equations. Hence, the procedureKlYa-JMMM:98  does not provide any gain with respect to former ones neither in analytical considerations nor in constructing fitting codes.

In the presented framework the results by Klik and Yao (including the analytical formulas for them missing in Ref. KlYa-JMMM:98, ) can be obtained immediately if to take the function $`\phi _1`$ in a step-wise form (23) and not to allow for the corrections caused by the finiteness of its derivative at $`x=0`$. In our terms this means to stop at set (45), i.e., “zero-derivative” solution, and not to go further. The emerging error is however, uncontrollable and not at all small. As an illustration, in Fig. 5 we show the result obtained with this model (dashed lines) for the cubic susceptibility $`\chi _{3\omega }^{(3)}`$ in a textured system where the particle common axis $`𝒏`$ is tilted under the angle $`\beta =\pi /3`$ to the probing field. One can see that deviations are substantial.

In Ref. RaSt-PRB:97,  we have proposed, although without rigorous justification, a formula for the cubic susceptibility of a random assembly

$$\chi _{3\omega }^{(3)}=\frac{1}{4}\chi _0^{(3)}\frac{(1+2S_2^2)(1i\omega \tau _1)}{45(1+i\omega \tau _1)(1+3i\omega \tau _1)},$$
(83)
that proved to be well adjusted for approximating the results of numerical calculations in all the temperature interval and also appeared to be good for fitting experimental data.SpFi-JMMM:01  Now we see that this very expression follows from Eqs. (77)–(79) if to expand the coefficients $`b_{i0}^{(3)}`$ up to the zeroth order with respect to $`\sigma ^1`$. This justifies Eq. (83) as a formula yielding a correct frequency dispersion of the cubic susceptibility of a random assembly at low temperatures. The cause of its applicability at high temperatures is the exponential dependence of $`\tau _1`$ on $`\sigma `$. Indeed, in the frequencies range $`\omega \tau _01`$, where we work, the condition $`\sigma 1`$ means $`\tau _1\tau _0`$, and all the dispersion factors in Eq. (83) drop out. This transforms expression (83) in a correct static susceptibility that is also a true result. To avoid any confusion we remark that Eq. (83) differs from formula for $`\chi _{3\omega }^{(3)}`$ given in Ref. RaSt-PRB:97  by the coefficient ($`1/45`$) due to the difference in definitions: in Ref. RaSt-PRB:97  it was included in $`\chi _0^{(3)}`$.

Applying the similar procedure to Eqs. (80)–(IV.3) we get the expression for the fifth-order susceptibility

$$\chi _{5\omega }^{(5)}=\frac{1}{16}\chi _0^{(5)}\frac{(2+12S_2^2+4S_2^3)}{945}\frac{1\frac{21}{8}i\omega \tau _1\frac{3}{4}\omega ^2\tau _1^2}{(1+i\omega \tau _1)(1+3i\omega \tau _1)(1+5i\omega \tau _1)},$$
(84)
that, following the example of the already tested Eq. (83), has high chances to be a good approximation for $`\chi _{5\omega }^{(5)}`$ in the whole temperature interval. As we have already ascertained in Ref. RaSt-PRB:97 , the best interpolation expression for the relaxation time in the susceptibility formulas is

$$\tau _1=\tau _D\frac{e^\sigma 1}{2\sigma }\left[\frac{\sigma }{1+\sigma }\sqrt{\frac{\sigma }{\pi }}+2^{\sigma 1}\right]^1,$$

proposed in Refs. CoCr-JMMM:94 ; CrCr-JAP:94 .

## VI Conclusions

A consistent procedure yielding the integral relaxation time and initial nonlinear susceptibilities for an assembly of non-interacting superparamagnetic particles is constructed in the low-to-moderate temperature range. Starting from the micromagnetic kinetic equation that describes intrinsic rotary diffusion of the particle magnetic moment, we obtain the results in an analytical form. They are presented as asymptotic series with respect to the dimensionless parameter $`\sigma `$ that is the uniaxial anisotropy barrier height scaled with temperature. High-order expansion terms are easily accessible that allows to achieve any desirable extent of accuracy. This is proven by comparison of the proposed approximation with the numerically-exact results. The susceptibilities contain angular dependencies that allow one to consider the particle assemblies with any extent of orientational texture—from perfectly aligned to random. The new formulas stand closer to reality than those for a two-level system and are to facilitate considerably both analytical and numerical calculations in the theory of superparamagnetic relaxation in single-domain particles.

###### Acknowledgements.

Partial financial support from the International Association for the Promotion of Cooperation with Scientists from the New Independent States of the Former Soviet Union (INTAS) under Grant No. 01–2341 and by Award No. PE–009–0 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF) is gratefully acknowledged. The package Maple V for PC used in the calculational work was obtained by Institute of Continuous Media Mechanics in the framework of the EuroMath Network and Services for the New Independent States - Phase II (EmNet/NIS/II) project funded by INTAS under grant No. IA-003.

## Appendix A Evaluation of the expansion coefficients for eigenfunctions $`\psi _0`$ and $`\psi _1`$

Both functions $`\psi _0`$ and $`\psi _1`$ are uniaxially symmetrical about the anisotropy axis $`𝒏`$ and can be expanded in the Legendre polynomial series, see Eq. (32):

$`\psi _0=\frac{1}{2}\underset{k=0}{}(2k+1)S_kP_k(x),k=0,2,4,\mathrm{},`$ (85)
$`\psi _1=\frac{1}{2}\underset{k=1}{}(2k+1)Q_kP_k(x),k=1,3,5,\mathrm{},`$
where in accordance with the parity properties of the eigenfunctions non-zero terms are

$`S_0=1,`$ $`S_k=(P_k(x)\left|\psi _0\right),k=2,4,\mathrm{},`$
$`Q_k=(P_k(x)\left|\psi _1\right),k=1,3,5,\mathrm{}.`$
Taking into account that $`\psi _1=\psi _0\phi _1`$, where $`\psi _0`$ in a finite form is given by Eq. (11), one arrives at the general formula

$$_k=(1/R)_0^1P_k(x)e^{\sigma x^2}𝑑x,$$
(87)
where $``$ is $`S_k`$ for even and is $`Q_k`$ for odd values of the index, and the function $`R(\sigma )`$ is defined by Eq. (11). In particular

$$Q_1=(1/R)_0^1xe^{\sigma x^2}𝑑x=\frac{1}{2}(e^\sigma 1)/\sigma R.$$
(88)
Using asymptotic expansion (12) for $`R`$, one gets

$$Q_1=1/G=1\frac{1}{2\sigma }\frac{1}{2\sigma ^2}\frac{5}{4\sigma ^3}\frac{37}{8\sigma ^4}\frac{353}{16\sigma ^5}\frac{4881}{32\sigma ^6}\frac{55205}{64\sigma ^7}\frac{854197}{128\sigma ^8}+\mathrm{}$$
(89)
Knowing $`Q_1`$, one can derive all the other moments $`Q_k`$ with the aid of the three-term recurrence relation obtained from Eq. (8) by setting there $`b_k=Q_k`$ and $`\lambda =0`$. The same relation can be used to find the equilibrium order parameters $`S_k`$. This is a head-to-tail procedure, where $`S_0=1`$ and $`S_2`$ is determined by the integral

$$S_2=(1/2R)_0^1\left(3x^21\right)e^{\sigma x^2}𝑑x.$$
(90)
Taking the latter by parts one gets

$$S_2=\frac{3}{4}\left[e^\sigma R\right]/\sigma R.$$

On comparison with Eq. (88), we find

$$S_2=\frac{3}{2}Q_1\frac{3}{4}(32\sigma )/\sigma ,$$

that upon substituting asymptotic series (89), transforms into

$$S_2=1\frac{3}{2\sigma }\frac{3}{4\sigma ^2}\frac{15}{8\sigma ^3}\frac{111}{16\sigma ^4}\frac{1059}{32\sigma ^5}\frac{12243}{64\sigma ^6}\frac{165615}{128\sigma ^7}\frac{2562591}{256\sigma ^8}+\mathrm{}$$
(91)
## Appendix B Evaluation of the correcting coefficients $`D_n`$ in a general case

Let us present the solution of Eq. (42) in the form

$$f_1^{(n)}=\psi _0g_1^{(n)}+u^{(n)},$$
(92)
where the functions $`g_1^{(n)}`$ are rendered by formulas (45) and are not corrected with respect to the derivative $`d\phi _1/dx`$. Substituting Eq. (92) in (42) and taking into account Eqs. (45), we get a recurrence sequence of equations for the corrections $`u^{(n)}`$:

$$\widehat{L}u^{(n)}=\widehat{V}u^{(n1)}+\widehat{𝑱}\psi _0\frac{(𝒆𝒉)^n}{n!}\widehat{𝑱}\phi _1.$$
(93)
With allowance for the fact that function $`\phi _1^{(0)}`$ depends only on $`x`$, Eq. (93) rewrites as

$$\widehat{L}u^{(n)}=\widehat{V}u^{(n1)}+\frac{d}{dx}\left[\psi _0(1x^2)\frac{(𝒆𝒉)^n}{n!}\frac{d\phi _1}{dx}\right].$$

Finally, making use of the relation

$$\frac{d\phi _1}{dx}=\frac{\lambda _1}{2\psi _0(1x^2)},$$
(94)
that follows from Eq. (18), we get

$$\widehat{L}u^{(n)}=\widehat{V}u^{(n1)}+\frac{\lambda _1}{2}\frac{d}{dx}\left[\frac{(𝒆𝒉)^n}{n!}\right].$$
(95)
In particular, at $`n=1`$ Eq. (95) takes the form

$$\widehat{L}u^{(1)}=\frac{\lambda _1}{2}\frac{d}{dx}(𝒆𝒉)$$
(96)
Equations (95) are solved sequentially beginning from Eq. (96) by expanding in a power series with respect to $`x`$. The right-hand sides of Eqs. (95) and (96) are proportional to an exponentially small parameter $`\lambda _1`$. Just due to that we did not take into account the corrections of the order $`u^{(n)}`$ when deriving Eqs. (45). However, the quantities

$$D_n=(\phi _1\left|u^{(n)}\right)n=2,4..,$$

have finite values. To show that, let us multiply Eq. (95) by $`\phi _1`$ and integrate. This yields

$$(\phi _1\left|\widehat{L}u^{(n)}\right)=(\phi _1\left|\widehat{V}u^{(n1)}\right)+\frac{\lambda _1}{2}(\phi _1\left|\frac{d}{dx}\left[\frac{(𝒆𝒉)^n}{n!}\right]\right).$$
(97)
In the left part we make use of the fact that $`\phi _1`$ is the left eigenfunction of the operator $`\widehat{L}`$, in the right part the integrals are taken by parts and yield

$`\lambda _1D_n=2{\displaystyle _0^1}(1x^2)u^{(n1)}\left[{\displaystyle \frac{d\phi _1}{dx}}\right]{\displaystyle \frac{d}{dx}}(𝒆𝒉)𝑑x`$ (98)
$`\lambda _1{\displaystyle _0^1}{\displaystyle \frac{d\phi _1}{dx}}{\displaystyle \frac{(𝒆𝒉)^n}{n!}}𝑑x.`$
Replacing the derivative $`d\phi _1/dx`$ in the first term of the right-hand side with the aid of Eq. (94), we arrive at the representation of the coefficient $`D_n`$ as

$$D_n=_0^1\frac{u^{(n1)}}{\psi _0}\frac{d}{dx}(𝒆𝒉)𝑑x_0^1\frac{d\phi _1}{dx}\frac{(𝒆𝒉)^n}{n!}𝑑x.$$
(99)
Since $`\psi _0\mathrm{exp}(\sigma x^2)`$, the first integral in Eq. (99) can be presented as an asymptotic series if the power expansion of the function $`u^{(n1)}`$ in the vicinity of $`x=0`$ is known. A closed form for the second integral can be found with the aid of the table given in Eq. (51), see section III.1.

As an example, we calculate the coefficient $`D_2`$. As from the addition theorem

$$(𝒆𝒉)=\mathrm{cos}\theta \mathrm{cos}\beta +\mathrm{sin}\theta \mathrm{sin}\beta \mathrm{cos}\phi ,$$

we seek the solution of Eq. (96) the sum

$$u^{(1)}=\mathrm{cos}\beta \underset{k}{}C_k^{(0)}x^k+\mathrm{sin}\beta e^{i\phi }(1x^2)^{\frac{1}{2}}\underset{k}{}C_k^{(1)}x^k.$$
(100)
Here the upper index of the $`C`$ coefficients corresponds to the azimuthal number $`m`$ of the spherical harmonic $`e^{im\phi }`$. Operator $`\widehat{L}`$ now includes the azimuthal coordinate and takes the form

$`\widehat{L}=(1x^2){\displaystyle \frac{d}{dx^2}}\left[2\sigma x(1x^2)+2x\right]{\displaystyle \frac{d}{dx}}`$
$`+\left[2\sigma (3x^21){\displaystyle \frac{m^2}{1x^2}}\right].`$
Substitution of expansion (100) in Eq. (96) leads to the set of equations

$$2\sigma (k+m+1)C_{k2}^{(m)}\left[k(k+1+2m+2\sigma )+m(m+1)+2\sigma \right]C_k^{(m)}+(k+1)(k+2)C_{k+2}^{(m)}=N_k^{(m)},$$
(101)
where $`m=0,1`$ and the numbers in the right-hand side are

$$N_k^{(0)}=\{\begin{array}{cc}1\hfill & \text{for}k=0,\hfill \\ 0\hfill & \text{for}k0,\hfill \end{array}N_k^{(1)}=\{\begin{array}{cc}1,\hfill & \text{for }k\text{ odd},\hfill \\ 0,\hfill & \text{for }k\text{ even}.\hfill \end{array}$$

In reality, one retains in expansion (100) only a finite number of terms so that Eqs. (101) could be easily solved analytically by any computer algebra solver. In terms of expansion (100) expression (99) at $`n=2`$ writes

$`D_2=\mathrm{cos}^2\beta {\displaystyle \underset{k=0}{}}C_{2k}^{(0)}{\displaystyle \frac{(2k1)!!}{2^k\sigma ^kG}}`$ (102)
$`{\displaystyle \frac{1}{2}}\mathrm{sin}^2\beta {\displaystyle \underset{k=1}{}}C_{2k1}^{(1)}{\displaystyle \frac{(2k1)!!}{2^k\sigma ^kG}}`$
$`{\displaystyle \frac{1}{6}}{\displaystyle \frac{2G3}{6G}}P_2(\mathrm{cos}\beta ).`$
Since the coefficients $`C`$ found from Eq. (101) are functions of $`\sigma `$, one has to perform in Eq. (102) asymptotic expansion. This gives finally

$`D_2={\displaystyle \frac{1}{4\sigma }}+{\displaystyle \frac{1}{4\sigma ^2}}+{\displaystyle \frac{5}{8\sigma ^3}}+{\displaystyle \frac{37}{16\sigma ^4}}+\mathrm{}`$ (103)
$`\mathrm{sin}^2\beta \left({\displaystyle \frac{1}{4}}+{\displaystyle \frac{1}{8\sigma }}+{\displaystyle \frac{1}{16\sigma ^2}}+{\displaystyle \frac{7}{64\sigma ^3}}+{\displaystyle \frac{19}{64\sigma ^4}}+\mathrm{}\right).`$
As it should be, at $`\beta =0`$ this formula reduces to Eq. (52) that was obtained for a one-dimensional case. We remark, however, that in a tilted situation ($`\beta 0`$) the coefficient $`D_2`$ acquires a contribution independent on $`\sigma `$ that assumes the leading role. This effect is clearly due to admixing of transverse modes to the set of eigenfunctions of the system, and it is just it that causes so a significant discrepancy between the “zero-derivative” approximation and the correct asymptotic expansion for $`\chi ^{(3)}`$ curves in Fig. 5. Evaluation of the coefficient $`D_4`$ is done according to the same scheme and requires taking into account a number of the perturbation terms that makes it rather cumbersome.

## Appendix C Evaluation of integrals

Before proceeding to the integrals (scalar products) in Eqs. (62)–(66) and (68), let us consider the “primitive” ones

$$X_n=((𝒆𝒉)^n\left|\psi _0\right),Y_n=((𝒆𝒉)^n\left|\psi _1\right).$$

The functions $`\psi _0`$ and $`\psi _1`$ are originally defined in terms of the angle $`\theta =\mathrm{arccos}(𝒆𝒏)`$. Thus, before performing integration one needs to transform both integrands to the same set of angles. Doing this with the aid of the addition theorem for Legendre polynomials, one finds

$`X_2=\frac{1}{3}[\mathrm{\hspace{0.17em}2}S_2P_2(\mathrm{cos}\beta )+1],X_4=\frac{1}{35}[\mathrm{\hspace{0.17em}8}S_4P_4(\mathrm{cos}\beta )+20S_2P_2(\mathrm{cos}\beta )+7],`$
$`X_6=\frac{1}{231}[\mathrm{\hspace{0.17em}16}S_6P_6(\mathrm{cos}\beta )+72S_4P_4(\mathrm{cos}\beta )+110S_2P_2(\mathrm{cos}\beta )+33],`$ (104)
and

$`Y_1=Q_1\mathrm{cos}\beta ,Y_3=\frac{1}{5}[\mathrm{\hspace{0.17em}2}Q_3P_3(\mathrm{cos}\beta )+3Q_1\mathrm{cos}\beta ],`$
$`Y_5=\frac{1}{63}[\mathrm{\hspace{0.17em}8}Q_5P_5(\mathrm{cos}\beta )+28Q_3P_3(\mathrm{cos}\beta )+27Q_1\mathrm{cos}\beta ],`$ (105)
where $`\mathrm{cos}\beta =(𝒏𝒉)`$ and the parameters $`S_k`$ and $`Q_k`$ are the expansion coefficients introduced by formulas (32).

Now using the expressions for functions $`f_0^{(n)}`$ and $`f_1^{(n)}`$ derived in Sec. III.1 one sees that the relevant integrals of Eqs. (63)–(68) are expressed in terms of $`X_k`$ and $`Y_k`$ as

$$((𝒆𝒉)\left|f_0^{(1)}\right)=X_2,((𝒆𝒉)\left|f_0^{(3)}\right)=\frac{1}{6}X_4\frac{1}{2}X_2^2,((𝒆𝒉)\left|f_0^{(5)}\right)=\frac{1}{120}X_6\frac{1}{8}X_4X_2+\frac{1}{4}X_2^3;$$
(106)
$$(\phi _1\left|f_0^{(1)}\right)=Y_1,(\phi _1\left|f_0^{(3)}\right)=\frac{1}{6}Y_3\frac{1}{2}X_2Y_1,(\phi _1\left|f_0^{(5)}\right)=\frac{1}{120}Y_5\frac{1}{12}Y_3X_2+\frac{1}{4}X_2^2Y_1\frac{1}{24}X_4Y_1;$$
(107)
$`(\phi _1\left|f_1^{(2)}\right)=\frac{1}{2}X_2Y_1^2+D_2,((𝒆𝒉)\left|f_1^{(2)}\right)=\frac{1}{2}Y_3Y_1X_2+D_2Y_1,`$
$`(\phi _1\left|f_1^{(4)}\right)=\frac{1}{24}X_4\frac{1}{3}Y_3Y_1+\frac{1}{2}X_2Y_1^2+D_4+D_2(\phi _1\left|f_1^{(2)}\right),`$
$`((𝒆𝒉)\left|f_1^{(4)}\right)=\frac{1}{24}Y_5\frac{1}{6}Y_3X_2\frac{1}{6}X_4Y_1+\frac{1}{2}X_2^2Y_1+D_4Y_1+D_2((𝒆𝒉)\left|f_1^{(2)}\right).`$ (108)


text:
# 1 Introduction

## 1 Introduction

Historically, variational principles have played a fundamental role in the evolution of mathematical models in classical physics, and many equations can be derived by using them. Here the relevant examples are Fermat’s principle in optics and Maupertuis’ principle in mechanics. The beginning of the analogy between geometrical optics and mechanics is usually attributed to Descartes (1637), but actually it can traced back to Ibn Al-Haitham Alhazen (0965-1037) . The analogy between the trajectory of material particles in potential fields and the path of light rays in media with continuously variable refractive index was formalized by Hamilton in 1833. This Hamiltonian analogy lead to the development of electron optics in 1920s, when Busch derived the focusing action and a lens-like action of the axially symmetric magnetic field using the methodology of geometrical optics. Around the same time Louis de Broglie associated his now famous wavelength to moving particles. Schrödinger extended the analogy by passing from geometrical optics to wave optics through his wave equation incorporating the de Broglie wavelength. This analogy played a fundamental role in the early development of quantum mechanics. The analogy, on the other hand, lead to the development of practical electron optics and one of the early inventions was the electron microscope by Ernst Ruska. A detailed account of Hamilton’s analogy is available in -.

Until very recently, it was possible to see this analogy only between the geometrical-optic and classical prescriptions of electron optics. The reasons being that, the quantum theories of charged-particle beam optics have been under development only for about a decade - with the very expected feature of wavelength-dependent effects, which have no analogue in the traditional descriptions of light beam optics. With the current development of the non-traditional prescriptions of Helmholtz optics  and the matrix formulation of Maxwell optics -, accompanied with wavelength-dependent effects, it is seen that the analogy between the two systems persists. The non-traditional prescription of Helmholtz optics is in close analogy with the quantum theory of charged-particle beam optics based on the Klein-Gordon equation. The matrix formulation of Maxwell optics is in close analogy with the quantum theory of charged-particle beam optics based on the Dirac equation. This analogy is summarized in the table of Hamiltonians. In this short note it is difficult to present the derivation of the various Hamiltonians which are available in the references. We shall briefly consider an outline of the quantum prescriptions and the non-traditional prescriptions respectively. A complete coverage to the new field of Quantum Aspects of Beam Physics (QABP), can be found in the proceedings of the series of meetings under the same name .

## 2 Quantum Formalism

The classical treatment of charged-particle beam optics has been extremely successful in the designing and working of numerous optical devices, from electron microscopes to very large particle accelerators. It is natural, however to look for a prescription based on the quantum theory, since any physical system is quantum mechanical at the fundamental level! Such a prescription is sure to explain the grand success of the classical theories and may also help get a deeper understanding and to lead to better designing of charged-particle beam devices.

The starting point to obtain a quantum prescription of charged particle beam optics is to build a theory based on the basic equations of quantum mechanics (Schrödinger, Klein-Gordon, Dirac) appropriate to the situation under study. In order to analyze the evolution of the beam parameters of the various individual beam optical elements (quadrupoles, bending magnets, $`\mathrm{}`$) along the optic axis of the system, the first step is to start with the basic time-dependent equations of quantum mechanics and then obtain an equation of the form

$$\mathrm{i}\mathrm{}\frac{}{s}\psi (x,y;s)=\widehat{}(x,y;s)\psi (x,y;s),$$
(1)
where $`(x,y;s)`$ constitute a curvilinear coordinate system, adapted to the geometry of the system. Eq. (1) is the basic equation in the quantum formalism, called as the beam-optical equation; $``$ and $`\psi `$ as the beam-optical Hamiltonian and the beam wavefunction respectively. The second step requires obtaining a relationship between any relevant observable $`\{O(s)\}`$ at the transverse-plane at $`s`$ and the observable $`\{O(s_{\mathrm{in}})\}`$ at the transverse plane at $`s_{\mathrm{in}}`$, where $`s_{\mathrm{in}}`$ is some input reference point. This is achieved by the integration of the beam-optical equation in (1)

$`\psi (x,y;s)`$ $`=`$ $`\widehat{U}(s,s_{\mathrm{in}})\psi (x,y;s_{\mathrm{in}}),`$ (2)
which gives the required transfer maps

$`O\left(s_{\mathrm{in}}\right)O\left(s\right)`$ $`=`$ $`\psi (x,y;s)\left|O\right|\psi (x,y;s),`$ (3)
$`=`$ $`\psi (x,y;s_{\mathrm{in}})\left|\widehat{U}^{}O\widehat{U}\right|\psi (x,y;s_{\mathrm{in}}).`$
The two-step algorithm stated above gives an over-simplified picture of the quantum formalism. There are several crucial points to be noted. The first step in the algorithm of obtaining the beam-optical equation is not to be treated as a mere transformation which eliminates $`t`$ in preference to a variable $`s`$ along the optic axis. A clever set of transforms are required which not only eliminate the variable $`t`$ in preference to $`s`$ but also give us the $`s`$-dependent equation which has a close physical and mathematical analogy with the original $`t`$-dependent equation of standard time-dependent quantum mechanics. The imposition of this stringent requirement on the construction of the beam-optical equation ensures the execution of the second-step of the algorithm. The beam-optical equation is such that all the required rich machinery of quantum mechanics becomes applicable to the computation of the transfer maps that characterize the optical system. This describes the essential scheme of obtaining the quantum formalism. The rest is mostly mathematical detail which is inbuilt in the powerful algebraic machinery of the algorithm, accompanied with some reasonable assumptions and approximations dictated by the physical considerations. The nature of these approximations can be best summarized in the optical terminology as a systematic procedure of expanding the beam optical Hamiltonian in a power series of $`|\widehat{𝝅}_{}/p_0|`$, where $`p_0`$ is the design (or average) momentum of beam particles moving predominantly along the direction of the optic axis and $`\widehat{𝝅}_{}`$ is the small transverse kinetic momentum. The leading order approximation along with $`|\widehat{𝝅}_{}/p_0|1`$, constitutes the paraxial or ideal behaviour and higher order terms in the expansion give rise to the nonlinear or aberrating behaviour. It is seen that the paraxial and aberrating behaviour get modified by the quantum contributions which are in powers of the de Broglie wavelength ($`\overline{\lambda }_0=\mathrm{}/p_0`$). The classical limit of the quantum formalism reproduces the well known Lie algebraic formalism of charged-particle beam optics .

## 3 Light Optics: Non-Traditional Prescriptions

The traditional scalar wave theory of optics (including aberrations to all orders) is based on the beam-optical Hamiltonian derived by using Fermat’s principle. This approach is purely geometrical and works adequately in the scalar regime. The other approach is based on the square-root of the Helmholtz operator, which is derived from the Maxwell equations . This approach works to all orders and the resulting expansion is no different from the one obtained using the geometrical approach of Fermat’s principle. As for the polarization: a systematic procedure for the passage from scalar to vector wave optics to handle paraxial beam propagation problems, completely taking into account the way in which the Maxwell equations couple the spatial variation and polarization of light waves, has been formulated by analyzing the basic Poincaré invariance of the system, and this procedure has been successfully used to clarify several issues in Maxwell optics -.

In the above approaches, the beam-optics and the polarization are studied separately, using very different machineries. The derivation of the Helmholtz equation from the Maxwell equations is an approximation as one neglects the spatial and temporal derivatives of the permittivity and permeability of the medium. Any prescription based on the Helmholtz equation is bound to be an approximation, irrespective of how good it may be in certain situations. It is very natural to look for a prescription based fully on the Maxwell equations, which is sure to provide a deeper understanding of beam-optics and light polarization in a unified manner.

The two-step algorithm used in the construction of the quantum theories of charged-particle beam optics is very much applicable in light optics! But there are some very significant conceptual differences to be borne in mind. When going beyond Fermat’s principle the whole of optics is completely governed by the Maxwell equations, and there are no other equations, unlike in quantum mechanics, where there are separate equations for, spin-$`1/2`$, spin-$`1`$, $`\mathrm{}`$.

Maxwell’s equations are linear (in time and space derivatives) but coupled in the fields. The decoupling leads to the Helmholtz equation which is quadratic in derivatives. In the specific context of beam optics, purely from a calculational point of view, the starting equations are the Helmholtz equation governing scalar optics and for a more accurate prescription one uses the full set of Maxwell equations, leading to vector optics. In the context of the two-step algorithm, the Helmholtz equation and the Maxwell equations in a matrix representation can be treated as the ‘basic’ equations, analogue of the basic equations of quantum mechanics. This works perfectly fine from a calculational point of view in the scheme of the algorithm we have.

Exploiting the similarity between the Helmholtz wave equation and the Klein-Gordon equation, the former is linearized using the Feshbach-Villars procedure used for the linearization of the Klein-Gordon equation. Then the Foldy-Wouthuysen iterative diagonalization technique is applied to obtain a Hamiltonian description for a system with varying refractive index. This technique is an alternative to the conventional method of series expansion of the radical. Besides reproducing all the traditional quasiparaxial terms, this method leads to additional terms, which are dependent on the wavelength, in the optical Hamiltonian. This is the non-traditional prescription of scalar optics.

The Maxwell equations are cast into an exact matrix form taking into account the spatial and temporal variations of the permittivity and permeability. The derived representation using $`8\times 8`$ matrices has a close algebraic analogy with the Dirac equation, enabling the use of the rich machinery of the Dirac electron theory. The beam optical Hamiltonian derived from this representation reproduces the Hamiltonians obtained in the traditional prescription along with wavelength-dependent matrix terms, which we have named as the polarization terms. These polarization terms are very similar to the spin terms in the Dirac electron theory and the spin-precession terms in the beam-optical version of the Thomas-BMT equation . The matrix formulation provides a unified treatment of beam optics and light polarization. Some well known results of light polarization are obtained as a paraxial limit of the matrix formulation -. The traditional beam optics is completely obtained from our approach in the limit of small wavelength, $`\text{}\lambda 0`$, which we call as the traditional limit of our formalisms. This is analogous to the classical limit obtained by taking $`\mathrm{}0`$, in the quantum prescriptions.

From the Hamiltonians in the Table we make the following observations: The classical/traditional Hamiltonians of particle/light optics are modified by wavelength-dependent contributions in the quantum/non-traditional prescriptions respectively. The algebraic forms of these modifications in each row is very similar. This should not come as a big surprise. The starting equations have one-to-one algebraic correspondence: Helmholtz $``$ Klein-Gordon; Matrix form of Maxwell $``$ Dirac equation. Lastly, the de Broglie wavelength, $`\text{}\lambda _0`$, and $`\lambda `$ have an analogous status, and the classical/traditional limit is obtained by taking $`\text{}\lambda _00`$ and $`\text{}\lambda 0`$ respectively. The parallel of the analogies between the two systems is sure to provide us with more insights.

## 4 Hamiltonians in Different Prescriptions

The following are the Hamiltonians, in the different prescriptions of light beam optics and charged-particle beam optics for magnetic systems. $`\widehat{H}_{0,p}`$ are the paraxial Hamiltonians, with lowest order wavelength-dependent contributions.

| Light Beam Optics | Charged-Particle Beam Optics |
| --- | --- |
| Fermat’s Principle $`=\left\{n^2\left(𝒓\right)𝒑_{}^2\right\}^{1/2}`$ | Maupertuis’ Principle $`=\left\{p_0^2𝝅_{}^2\right\}^{1/2}qA_z`$ |
| Non-Traditional Helmholtz $`\widehat{H}_{0,p}=n\left(𝒓\right)+\frac{1}{2n_0}\widehat{𝒑}_{}^2\frac{\mathrm{i}\text{}\lambda }{16n_0^3}[\widehat{𝒑}_{}^2,\frac{}{z}n\left(𝒓\right)]`$ | Klein-Gordon Formalism $`\widehat{H}_{0,p}=p_0qA_z+\frac{1}{2p_0}\widehat{𝝅}_{}^2+\frac{\mathrm{i}\mathrm{}}{16p_0^4}[\widehat{𝝅}_{}^2,\frac{}{z}\widehat{𝝅}_{}^2]`$ |
| Maxwell, Matrix $`\widehat{H}_{0,p}=n\left(𝒓\right)+\frac{1}{2n_0}\widehat{𝒑}_{}^2\mathrm{i}\text{}\lambda \beta 𝚺𝒖\text{}+\frac{1}{2n_0}\text{}\lambda ^2w^2\beta `$ | Dirac Formalism $`\widehat{H}_{0,p}=p_0qA_z+\frac{1}{2p_0}\widehat{𝝅}_{}^2\frac{\mathrm{}}{2p_0}\left\{\mu \gamma 𝚺_{}𝑩_{}+\left(q+\mu \right)\mathrm{\Sigma }_zB_z\right\}+\mathrm{i}\frac{\mathrm{}}{m_0c}ϵB_z`$ |

Notation

| $`\mathrm{Refractive}\mathrm{Index},n\left(𝒓\right)=c\sqrt{ϵ\left(𝒓\right)\mu \left(𝒓\right)}\mathrm{Resistance},h\left(𝒓\right)=\sqrt{\mu \left(𝒓\right)/ϵ\left(𝒓\right)}𝒖\left(𝒓\right)=\frac{1}{2n\left(𝒓\right)}\mathbf{}n\left(𝒓\right)𝒘\left(𝒓\right)=\frac{1}{2h\left(𝒓\right)}\mathbf{}h\left(𝒓\right)`$ $`𝚺`$ and $`\beta `$ are the Dirac matrices. | $`\widehat{𝝅}_{}=\widehat{𝒑}_{}q𝑨_{}\mu _a\mathrm{anomalous}\mathrm{magnetic}\mathrm{moment}.ϵ_a\mathrm{anomalous}\mathrm{electric}\mathrm{moment}.\mu =2m_0\mu _a/\mathrm{},ϵ=2m_0ϵ_a/\mathrm{}\gamma =E/m_0c^2`$ |
| --- | --- |