text:
# Successive Minima and Lattice Points

## 1. Introduction

In 1896 Hermann Minkowski’s fundamental and guiding book “Geometrie der Zahlen” \[Min96\] was published, which may be considered as the first systematic study on relations between convex geometry, Diophantine approximation, and the theory of quadratic forms (cf. Gruber \[Gru93\]). One of the basic problems in geometry of numbers is to decide whether a given set in the $`d`$-dimensional Euclidean space $`^d`$ contains a non-trivial lattice point of a $`d`$-dimensional lattice $`\mathrm{\Lambda }^d`$. With respect to the class $`𝒦_0^d`$ of all 0-symmetric convex bodies in $`^d`$ with non-empty interior and the volume $`\mathrm{vol}()`$$`d`$-dimensional Lebesgue measure – Minkowski settled this problem:

(1.1) If $`\mathrm{vol}(K)2^ddet\mathrm{\Lambda }`$ then $`K`$ contains a non-zero lattice point of $`\mathrm{\Lambda }`$.
Here $`det\mathrm{\Lambda }`$ denotes the determinant of the lattice $`\mathrm{\Lambda }`$ and the space of all lattices $`\mathrm{\Lambda }^d`$ with $`det\mathrm{\Lambda }0`$ is denoted by $`^d`$.

Minkowski assessed his result as “ein Satz, der nach meinem Dafürhalten zu den fruchtbarsten in der Zahlenlehre zu rechnen ist” (\[Min96\], p. 75) and indeed this theorem has many applications (cf. \[EGH89\], sec. 3.3). Minkowski proved even a stronger result, for which we have to introduce his “kleinstes System von unabhängig gerichteten Strahlendistanzen im Zahlengitter” (\[Min96\], p. 178).

###### Definition 1.1.

Let $`K𝒦_0^d`$ and $`\mathrm{\Lambda }^d`$. For $`1id`$

$$\begin{array}{cc}\hfill \lambda _i(K,\mathrm{\Lambda })=\mathrm{min}\{\lambda _0:& \lambda K\text{contains }i\text{ linearly independent}\hfill \\ & \text{lattice points of }\mathrm{\Lambda }\}\hfill \end{array}$$

is called the $`i`$-th successive minimum of $`K`$ with respect to $`\mathrm{\Lambda }`$.

Obviously, we have $`\lambda _1(K,\mathrm{\Lambda })\lambda _2(K,\mathrm{\Lambda })\mathrm{}\lambda _d(K,\mathrm{\Lambda })`$ and the first successive minimum $`\lambda _1(K,\mathrm{\Lambda })`$ is the smallest dilation factor such that $`\lambda _1(K,\mathrm{\Lambda })K`$ contains a non-zero lattice point. With this notation Minkowski’s first theorem on successive minima reads (cf. \[Min96\], pp. 75)

###### Theorem 1.2 (Minkowski).

Let $`K`$ in $`𝒦_0^d`$ and $`\mathrm{\Lambda }^d`$. Then

$$\lambda _1(K,\mathrm{\Lambda })^d\mathrm{vol}(K)2^ddet\mathrm{\Lambda }.$$

So $`\mathrm{vol}(K)2^ddet\mathrm{\Lambda }`$ implies $`\lambda _1(K,\mathrm{\Lambda })1`$, and we get (1.1). Minkowski’s second theorem on successive minima is a deep improvement of the first one and says (cf. \[Min96\], pp. 199)

###### Theorem 1.3 (Minkowski).

Let $`K𝒦_0^d`$ and $`\mathrm{\Lambda }^d`$. Then

$$\lambda _1(K,\mathrm{\Lambda })\lambda _2(K,\mathrm{\Lambda })\mathrm{}\lambda _d(K,\mathrm{\Lambda })\mathrm{vol}(K)2^ddet\mathrm{\Lambda }.$$

This inequality is best possible. For instance, with respect to the integral lattice $`^d`$, each box with axes parallel to the coordinate axes gives equality. Although Theorem 1.3 has not so many applications as the first theorem on successive minima, it shows a beautiful relation between the volume of $`K`$ and the expansion of $`K`$ with respect to independent lattice directions of a lattice. The importance of Theorem 1.3 is also reflected in the number of different proofs, see e.g. Bambah, Woods & Zassenhaus \[BWZ65\], Cassels \[Cas59\], Danicic \[Dan69\], Davenport \[Dav39\], Estermann \[Est46\], Siegel \[Sie89\] and Weyl \[Wey42\].

In \[BHW93\] it was conjectured that an inequality analogue to Theorem 1.3 holds for the lattice point enumerator $`\mathrm{\#}(K\mathrm{\Lambda })`$. More precisely,

###### Conjecture 1.4.

Let $`K𝒦_0^d`$ and $`\mathrm{\Lambda }^d`$. Then

(1.2)
$$\mathrm{\#}(K\mathrm{\Lambda })\underset{i=1}{\overset{d}{}}\frac{2}{\lambda _i(K,\mathrm{\Lambda })}+1.$$

Here $`x`$ denotes the smallest integer not less than $`x`$. An analogous statement to the first theorem of Minkowski on successive minima was already shown in \[BHW93\], namely

(1.3)
$$\mathrm{\#}(K\mathrm{\Lambda })\frac{2}{\lambda _1(K,\mathrm{\Lambda })}+1^d.$$

It seems to be worth mentioning that if Conjecture 1.2 were true then we could write by the definition of the Riemann integral

$$\frac{\mathrm{vol}(K)}{det\mathrm{\Lambda }}=\underset{r0}{lim}r^d\mathrm{\#}(Kr\mathrm{\Lambda })\underset{r0}{lim}\underset{i=1}{\overset{d}{}}r\frac{2}{\lambda _i(K,r\mathrm{\Lambda })}+1=\underset{i=1}{\overset{d}{}}\frac{2}{\lambda _i(K,\mathrm{\Lambda })}.$$

Thus Conjecture 1.2 implies Minkowski’s second theorem on successive minima (Theorem 1.3). In \[BHW93\] the validity of the conjecture was proven in the case $`d=2`$. Moreover, it was shown that an upper bound of this type exists, if in the above product $`\frac{2}{\lambda _i(K,\mathrm{\Lambda })}`$ is replaced by $`\frac{2i}{\lambda _i(K,\mathrm{\Lambda })}`$. So, roughly speaking, (1.2) holds up to a factor $`d!`$. Here we shall improve this bound.

###### Theorem 1.5.

Let $`d2`$, $`K𝒦_0^d`$ and $`\mathrm{\Lambda }^d`$. Then

$$\mathrm{\#}(K\mathrm{\Lambda })<2^{d1}\underset{i=1}{\overset{d}{}}\frac{2}{\lambda _i(K,\mathrm{\Lambda })}+1.$$

The proof of this theorem will be given is the next section. Minkowski’s original proof (\[Min96\], 199-218) of his second theorem on successive minima was sometimes criticised as lengthy and obscure (cf. \[Dav77\], p.91). One reason might be that in the scope of the proof he also proves many basic facts about the volume of a convex body, like the computation of the volume through successive integrations, etc., which cloud a little bit the simple and nice geometrical ideas of his proof. Based on these ideas we present a short proof of Theorem 1.3 in the last section.

For more information on lattices, successive minima and their role in the geometry of numbers we refer to the books of Erdös, Gruber and Hammer \[EGH89\], Gruber and Lekkerkerker \[GL87\] and the survey of Gruber \[Gru93\]. For an elementary introduction to the geometry of numbers see \[OLD00\].

## 2. Proof of Theorem 1.5

Before giving the proof we list some basic facts on lattices, for which we refer to \[GL87\]. Every lattice $`\mathrm{\Lambda }^d`$ can be written as $`\mathrm{\Lambda }=A^d`$, where $`A`$ is a non-singular $`(d\times d)`$-matrix, i.e., $`A\mathrm{GL}(d,)`$. In particular we have $`\lambda _i(K,\mathrm{\Lambda })=\lambda _i(A^1K,^d)`$ and $`\mathrm{\#}(K\mathrm{\Lambda })=\mathrm{\#}(A^1K^d)`$. A lattice $`\stackrel{~}{\mathrm{\Lambda }}^d`$ is called a sublattice of $`\mathrm{\Lambda }^d`$ if $`\stackrel{~}{\mathrm{\Lambda }}\mathrm{\Lambda }`$. For $`a,\overline{a}\mathrm{\Lambda }`$ and a sublattice $`\stackrel{~}{\mathrm{\Lambda }}\mathrm{\Lambda }`$ we write

$$a\overline{a}mod\stackrel{~}{\mathrm{\Lambda }}(a\overline{a})\stackrel{~}{\mathrm{\Lambda }}.$$

In words, $`a,\overline{a}`$ belong to the same residue class (coset) of $`\mathrm{\Lambda }`$ with respect to $`\stackrel{~}{\mathrm{\Lambda }}`$. We note that there are precisely $`det\stackrel{~}{\mathrm{\Lambda }}/det\mathrm{\Lambda }`$ different residue classes of $`\mathrm{\Lambda }`$ with respect to $`\stackrel{~}{\mathrm{\Lambda }}`$. For every set of $`d`$-linearly independent lattice points $`a^1,\mathrm{},a^d`$ of a lattice $`\mathrm{\Lambda }`$ there exists a basis $`b^1,\mathrm{},b^d`$ of $`\mathrm{\Lambda }`$ such that $`\mathrm{lin}\{a^1,\mathrm{},a^i\}=\mathrm{lin}\{b^1,\mathrm{},b^i\}`$, where $`\mathrm{lin}`$ denotes the linear hull. In particular, given $`d`$ linearly independent lattice vectors $`z^i^d`$, $`1id`$, with $`z^i\lambda _i(K,^d)K`$ then there exists an unimodular matrix $`U`$, i.e., $`U\mathrm{GL}(d,)^{d\times d}`$, such that

(2.1)
$$Uz^i\left(\lambda _i(UK,^d)UK\right)\mathrm{lin}\{e^1,\mathrm{},e^i\},1id,$$

where $`e^i^d`$ denotes the $`i`$-th unit vector. Furthermore we note that for $`d`$ linearly independent lattice points $`a^1,\mathrm{},a^d`$ of a lattice $`\mathrm{\Lambda }^d`$ satisfying $`a^i\lambda _i(K,\mathrm{\Lambda })K`$, the definition of the successive minima implies

(2.2)
$$\mathrm{int}\left(\lambda _i(K,\mathrm{\Lambda })K\right)\mathrm{\Lambda }\mathrm{lin}\{0,a^1,\mathrm{},a^{i1}\}\mathrm{\Lambda },1id,$$

where $`\mathrm{int}`$ denotes the interior.

For the proof of Theorem 1.5 we need the following simple lemma.

###### Lemma 2.1.

Let $`K𝒦_0^d`$, $`\mathrm{\Lambda }^d`$ and let $`\stackrel{~}{\mathrm{\Lambda }}`$ be a sublattice of $`\mathrm{\Lambda }`$. Then

$$\mathrm{\#}\left(K\mathrm{\Lambda }\right)\frac{det\stackrel{~}{\mathrm{\Lambda }}}{det\mathrm{\Lambda }}\mathrm{\#}\left(2K\stackrel{~}{\mathrm{\Lambda }}\right).$$

###### Proof.

Let $`m=\mathrm{\#}(2K\stackrel{~}{\mathrm{\Lambda }})`$ and suppose there exist at least $`m+1`$ different lattice points $`a^1,\mathrm{},a^{m+1}K\mathrm{\Lambda }`$ such that $`a^ia^1mod\stackrel{~}{\mathrm{\Lambda }}`$, $`1im+1`$. Then we have

$$a^ia^1(KK)\stackrel{~}{\mathrm{\Lambda }}=2K\stackrel{~}{\mathrm{\Lambda }},1im+1,$$

which contradicts the assumption $`\mathrm{\#}(2K\stackrel{~}{\mathrm{\Lambda }})=m`$. Thus we have shown that every residue class of $`\mathrm{\Lambda }`$ with respect to $`\stackrel{~}{\mathrm{\Lambda }}`$ does not contain more than $`m`$ points of $`K\mathrm{\Lambda }`$. Since there are precisely $`det\stackrel{~}{\mathrm{\Lambda }}/det\mathrm{\Lambda }`$ different residue classes, we get the desired bound. ∎

We remark that inequality (1.3) is a simple consequence of this lemma. To see this we set $`n_1=2/\lambda _1(K,\mathrm{\Lambda })+1`$ and $`\stackrel{~}{\mathrm{\Lambda }}=n_1\mathrm{\Lambda }`$. Next we observe that $`a2K\stackrel{~}{\mathrm{\Lambda }}`$ implies that $`\frac{1}{n_1}a\mathrm{int}(\lambda _1(K,\mathrm{\Lambda })K)\mathrm{\Lambda }`$ and from (2.2) we conclude $`\mathrm{\#}(2K\stackrel{~}{\mathrm{\Lambda }})=1`$. Thus Lemma 2.1 gives

$$\mathrm{\#}(K\mathrm{\Lambda })\frac{det\stackrel{~}{\mathrm{\Lambda }}}{det\mathrm{\Lambda }}=(n_1)^d=\frac{2}{\lambda _1(K,\mathrm{\Lambda })}+1^d.$$

Next we come to the proof of Theorem 1.5.

###### Proof of Theorem 1.5.

W.l.o.g. let $`\mathrm{\Lambda }=^d`$ and we may assume that (cf. (2.1) and (2.2))

(2.3)
$$\mathrm{int}\left(\lambda _i(K,^d)K\right)^d\mathrm{lin}\{0,e^1,\mathrm{},e^{i1}\}^d,1id.$$

For abbreviation we set $`q_i=\frac{2}{\lambda _i(K,\mathrm{\Lambda })}+1`$, $`1id`$, and first we determine $`d`$ numbers $`n_i`$ such that

(2.4)
$$\begin{array}{c}\hfill n_d=q_d,q_in_i<2q_i,\text{ and }n_{i+1}\text{ divides }n_i,\mathrm{\hspace{0.17em}1}id1.\end{array}$$

Suppose we have already found $`n_d,\mathrm{},n_{k+1}`$ with these properties. In order to determine $`n_k`$ we distinguish two cases. If $`n_{k+1}q_k`$ we set $`n_k=n_{k+1}`$. Since $`q_kq_{k+1}`$ we obtain $`q_kn_k=n_{k+1}<2q_{k+1}2q_k`$. Otherwise, if $`n_{k+1}<q_k`$ let $`q_k=mn_{k+1}+r`$ with $`m`$, $`m1`$, and $`0r<n_{k+1}`$. In this case we set $`n_k=q_k+n_{k+1}r`$ and obviously, $`n_k`$ meets the requirements of (2.4).

Now let $`\stackrel{~}{\mathrm{\Lambda }}^d`$ be the lattice generated by the vectors $`n_1e^1,n_2e^2,\mathrm{},n_de^d`$. Then we have $`det\stackrel{~}{\mathrm{\Lambda }}/det\mathrm{\Lambda }=n_1n_2\mathrm{}n_d`$ and together with the upper bounds on the the numbers $`n_i`$, Lemma 2.1 gives

(2.5)
$$\mathrm{\#}\left(K\mathrm{\Lambda }\right)\mathrm{\#}\left(2K\stackrel{~}{\mathrm{\Lambda }}\right)\underset{i=1}{\overset{d}{}}n_i<\mathrm{\#}\left(2K\stackrel{~}{\mathrm{\Lambda }}\right)\mathrm{\hspace{0.17em}2}^{d1}\underset{i=1}{\overset{d}{}}\frac{2}{\lambda _i(K,\mathrm{\Lambda })}+1.$$

Hence, in order to verify the theorem, it suffices to show $`2K\stackrel{~}{\mathrm{\Lambda }}=\{0\}`$. Suppose there exists a $`g2K\stackrel{~}{\mathrm{\Lambda }}\{0\}`$ and let $`k`$ be the largest index of a non-zero coordinate of $`g`$, i.e., $`g_k0`$ and $`g_{k+1}=\mathrm{}=g_d=0`$. Then we may write

$$g=z_1(n_1e^1)+z_2(n_2e^2)+\mathrm{}+z_k(n_ke^k)2K$$

for some $`z_i`$. Since $`n_k`$ is a divisor of $`n_1,\mathrm{},n_{k1}`$ and since $`2/n_k<\lambda _k(K,^d)`$ (cf. (2.4)) we obtain

$$\frac{1}{n_k}g\left(\frac{2}{n_k}K\right)^d\mathrm{int}(\lambda _k(K,^d)K)^d.$$

However, since $`g_k0`$ this relation violates (2.3). Thus we have $`2K\stackrel{~}{\mathrm{\Lambda }}=\{0\}`$ and the theorem is proven. ∎

## 3. Proof of Theorem 1.3

Minkowski’s proof of his second theorem on successive minima can be found in his book “Geometrie der Zahlen” (\[Min96\], 199–219) and for an English translation we refer to \[Han64\], 570–603.

###### Proof of Theorem 1.3 (following Minkowski).

Again w.l.o.g. we may assume that $`\mathrm{\Lambda }=^d`$. For convenience we write $`\lambda _i=\lambda _i(K,^d)`$ and set $`K_i=\frac{\lambda _i}{2}K`$. Furthermore, we assume that $`z^1,\mathrm{},z^d`$ are $`d`$ linearly independent lattice points with $`z^i\lambda _iK^d`$ and $`\mathrm{lin}\{z^1,\mathrm{},z^i\}=\mathrm{lin}\{e^1,\mathrm{},e^i\}`$, $`1id`$, (cf. (2.1)). For short, we denote the linear space $`\mathrm{lin}\{e^1,\mathrm{},e^i\}`$ by $`L_i`$.

For an integer $`q`$ let $`M_q^d=\{z^d:|z_i|q,\mathrm{\hspace{0.17em}1}id\}`$ and for $`1jd1`$ let $`M_q^j=M_q^dL_j`$. Since $`K`$ is a bounded set there exists a constant $`\gamma `$, only depending on $`K`$, such that

(3.1)
$$\mathrm{vol}(M_q^d+K_d)(2q+\gamma )^d.$$

By the definition of $`\lambda _1`$ we have $`(z+\mathrm{int}(K_1))(\overline{z}+\mathrm{int}(K_1))=\mathrm{}`$ for two different lattice point $`z,\overline{z}^d`$, because otherwise we would get the contradiction $`z\overline{z}(\mathrm{int}(K_1)\mathrm{int}(K_1))^d=\mathrm{int}(K_1K_1)^d=\mathrm{int}(\lambda _1K)^d=\{0\}`$. Thus we have

(3.2)
$$\mathrm{vol}(M_q^d+K_1)=(2q+1)^d\mathrm{vol}(K_1)=(2q+1)^d\left(\frac{\lambda _1}{2}\right)^d\mathrm{vol}(K).$$

In the following we shall show that for $`1id1`$

(3.3)
$$\mathrm{vol}(M_q^d+K_{i+1})\left(\frac{\lambda _{i+1}}{\lambda _i}\right)^{di}\mathrm{vol}(M_q^d+K_i).$$

To this end we may assume $`\lambda _{i+1}>\lambda _i`$ and let $`z,\overline{z}^d`$, which differ in the last $`di`$ coordinates, i.e., $`(z_{i+1},\mathrm{},z_d)(\overline{z}_{i+1},\mathrm{},\overline{z}_d)`$. Then

(3.4)
$$\left[z+\mathrm{int}(K_{i+1})\right]\left[\overline{z}+\mathrm{int}(K_{i+1})\right]=\mathrm{}.$$

Otherwise the $`i+1`$ linearly independent lattice points $`z\overline{z},z^1,\mathrm{},z^i`$ belong to the interior of $`\lambda _{i+1}K`$ which contradicts the minimality of $`\lambda _{i+1}`$. Hence we obtain from (3.4)

$$\begin{array}{cc}\hfill \mathrm{vol}\left(M_q^d+K_{i+1}\right)& =(2q+1)^{di}\mathrm{vol}\left(M_q^i+K_{i+1}\right),\hfill \\ \hfill \mathrm{vol}\left(M_q^d+K_i\right)& =(2q+1)^{di}\mathrm{vol}\left(M_q^i+K_i\right).\hfill \end{array}$$

and in order to verify (3.3) it suffices to show

(3.5)
$$\mathrm{vol}\left(M_q^i+K_{i+1}\right)\left(\frac{\lambda _{i+1}}{\lambda _i}\right)^{di}\mathrm{vol}(M_q^i+K_i).$$

Let $`f_1,f_2:^d^d`$ be the linear maps given by

$`f_1(x)`$ $`=`$ $`({\displaystyle \frac{\lambda _{i+1}}{\lambda _i}}x_1,\mathrm{},{\displaystyle \frac{\lambda _{i+1}}{\lambda _i}}x_i,x_{i+1},\mathrm{},x_d)^{},`$
$`f_2(x)`$ $`=`$ $`(x_1,\mathrm{},x_i,{\displaystyle \frac{\lambda _{i+1}}{\lambda _i}}x_{i+1},\mathrm{},{\displaystyle \frac{\lambda _{i+1}}{\lambda _i}}x_d)^{}.`$
Since $`M_q^i+K_{i+1}=f_2(M_q^i+f_1(K_i))`$ we get

$$\mathrm{vol}\left(M_q^i+K_{i+1}\right)=\left(\frac{\lambda _{i+1}}{\lambda _i}\right)^{di}\mathrm{vol}(M_q^i+f_1(K_i))$$

and for the proof of (3.5) we have to show

(3.6)
$$\mathrm{vol}\left(M_q^i+f_1(K_i)\right)\mathrm{vol}\left(M_q^i+K_i\right).$$

To this end let $`L_i^{}`$ be the $`(di)`$-dimensional orthogonal complement of $`L_i`$. Then it is easy to see that for every $`xL_i^{}`$ there exists a $`t(x)L_i`$ with $`K_i(x+L_i)(f_1(K_i)(x+L_i))+t(x)`$ and so

$$\left(M_q^i+K_i\right)\left(x+L_i\right)\left[\left(M_q^i+f_1(K_i)\right)\left(x+L_i\right)\right]+t(x).$$

Thus we get

$`\mathrm{vol}(M_q^i+K_i)`$ $`=`$ $`{\displaystyle _{xL_i^{}}}\mathrm{vol}_i\left((M_q^i+K_i)(x+L_i)\right)dx`$
$``$ $`{\displaystyle _{xL_i^{}}}\mathrm{vol}_i\left((M_q^i+f_1(K_i))(x+L_i)\right)dx`$
$`=`$ $`\mathrm{vol}(M_q^i+f_1(K_i)),`$
where $`\mathrm{vol}_i()`$ denotes the $`i`$-dimensional volume. This shows (3.6) and so we have verified (3.3). Finally, it follows from (3.1), (3.2) and (3.3)

$$\begin{array}{cc}\hfill (2q+\gamma )^d& \mathrm{vol}\left(M_q^d+K_d\right)\left(\frac{\lambda _d}{\lambda _{d1}}\right)\mathrm{vol}\left(M_q^d+K_{d1}\right)\hfill \\ & \left(\frac{\lambda _d}{\lambda _{d1}}\right)\left(\frac{\lambda _{d1}}{\lambda _{d2}}\right)^2\mathrm{vol}\left(M_q^d+K_{d2}\right)\mathrm{}\mathrm{}\hfill \\ & \left(\frac{\lambda _d}{\lambda _{d1}}\right)\left(\frac{\lambda _{d1}}{\lambda _{d2}}\right)^2\mathrm{}\left(\frac{\lambda _2}{\lambda _1}\right)^{d1}\mathrm{vol}\left(M_q^d+K_1\right)\hfill \\ & =\lambda _d\mathrm{}\lambda _1\frac{\mathrm{vol}(K)}{2^d}(2q+1)^d\hfill \end{array}$$

and so

$$\lambda _1\mathrm{}\lambda _d\mathrm{vol}(K)2^d\left(\frac{2q+\gamma }{2q+1}\right)^d.$$

Since this holds for all $`q`$ the theorem is proven. ∎


format: text

text:
# Magnetic dot arrays modeling via the system of the radial basis function networks

## 1 Introduction

In the recent years there is a remarkable progress in the technology of the nanofabrication of well defined magnetic materials. The material nanoscience based on the epitaxial and lithographic techniques  allows the fabrication of the regular arrays of the magnetic particles-dots of well controlled and interesting shape , lattice geometry and composition. The increasing technological flexibility calls for further physical ideas, which should be incorporated into design of the artificial nanoscale magnetic systems.

The uniformity of the polarization is the general basic aspect discussed in the connection with the small magnetic particles. The concept of uniformly polarized particle is justified only for the particles of an intermediate size . In the theory  the magnetostatic coupling was derived for the homogeneously polarized and saturated cylindrical dots on a rectangular lattice. More restrictive are conditions of the simulation , where each dot of array is substituted by a single dipolar moment. This approximation can be used only for monodomain dots separated by a sufficiently large distances. When a dot array is represented by a system of the interacting dipoles, the search for the ground state configuration leads to the formulation typical for the classical dipolar lattices . The violation of the intradot homogeneity stems from the competition between the magnetostatic, anisotropy and exchange energy terms. The analytical model of the dot array going towards the non-uniformity was proposed in . In this model the interactions of dots were described by the quadrupolar terms.

The problem of the calculation of the magnetization field of a dot array can be in principal formulated in the terms of classical micromagnetic theory . Due to complexity of the problem, the important role in its treatment will play the numerical simulations. They require the implementation of the sufficiently dense discretization within the each ferromagnetic dot. As usual, the magnetic part of the system can be subdivided into interacting dipoles or grains , small ferromagnetic cubes  or finite elements . Then the optimum spacing of the mesh nodes is determined by a minimum magnetic length scale (exchange, wall) of the system. For the majority of ferromagnetic materials, the comprehensive micromagnetic description is attained when the size of discretization elements decreases into the nanometer regime . Thus, the simulation of a single ultrathin dot of the micrometer size requires about $`10^6`$ nodes, although, the qualitative simulations can be realized even for $`10^210^4`$ nodes . From this we can conclude that both detailed and truncated micromagnetic description of many-dot array represents rather demanding computational task. We summarize, that principal difficulties of micromagnetic dot array analysis come from: (i) the interplay of the phenomena on the intradot (exchange, domain wall) and external geometric length scales; (ii) complexity of the magnetic structure of the non-uniformly polarized dots; (iii) long-range magnetostatic interdot interactions.

To make the problem of the magnetic ordering of dot arrays tractable by a moderate computer facilities, we have developed method, which works on a much coarser mesh than usual discretization schemes allow (except the adaptive and multigrid methods). Its general idea is the simultaneous simulation of the intradot - micromagnetic and multidot scales. This idea was strongly inspired by the multiscale approach . At the present stage of the project, the multidot part of simulation has been developed separately and the behaviour of the small-intradot scales has been treated only phenomenologically. The approach allows a remarkable increase of the simulation speed, indeed, the price payed to the scale separation is the appearance of additional parameters. The completing of the project needs support of the algorithms of the parameter estimation developed on the basis of the standard micromagnetic simulations. Let us to note that similar problems were solved in a cellular automaton version of the molecular dynamics .

In this paper we presented the results of the simulation of magnetic properties of quasi-two-dimensional cluster of ultra-thin magnetic dots on the square lattice. The phenomenological aspect of our model is a variable intradot inhomogeneity. The formalism we developed for this aim has been adopted from the models of the neural networks.

From the point of view of the information theory, the neural networks are continuous, unique mappings constructed from the system of known activation functions. The synaptic weights of these activation functions are adjusted by the training process. The standard problem, which can be effectively solved by the neural networks is the association of the input patterns (in our case inputs are effective magnetic moments) with the desired outputs (magnetic self-energy of dot). There are many applications, where neural networks can be implemented. They allow interpolation of the data generated by the simulations or experiment. The example of the physical application is , where neural network was used to fit a complicated analytic potential to the set of ab initio data. In , the Hopfield type of the interaction matrix was suggested to simulate the dynamics of the complex protein molecule. The specific magnetic application represents the solution of the magnetic inverse problem .

For the purpose to model the variable magnetic intradot inhomogeneity we adopted the theory of the radial basis function networks (RBFN) . The RBFN variant of the neural network was chosen, because its ingredient is a straightforward and explicit estimation of the synaptic weights, which allows more transparent analysis of the physical symmetries.

The aim of the paper is to make general presentation of the model and present some numerical results. The paper is organized as follows: In Section 2 the model of the dot array energy functional is introduced. This Section consists of two parts: in Subsection 2.1 we introduced the general phenomenological concept of the nonuniform magnetization, which utilizes self-energy interpolation by means of RBFN approach. In Subsection 2.2 the interaction between dots is introduced. In Section 3 our method is applied to the ultrathin square dots, where tendency of the formation of the vortex intradot phase prevails. Section 4 provides some details about the implementation of the simulated annealing algorithm to the problem of the total energy minimization. Finally, in Section 5 we bring examples of the numerical simulations.

## 2 Model

### 2.1 Intradot self-energy

The microstate of the system of $`N`$ magnetic dots is described by $`N_\mathrm{c}byN`$ effective magnetic moments

$$𝐦_{in},i=1,2,\mathrm{},N,n=1,2,\mathrm{},N_\mathrm{c},$$
(1)
which are associated with the magnetization field of dot $`𝐌(x,y)`$ via the volume averages

$$𝐦_{in}=m_{\mathrm{sat}}^1_{(x,y)\mathrm{}_{in}}dxdy𝐌(x,y).$$
(2)
The formal integration is performed here over $`n`$th volume element of $`i`$th dot labeled as $`\mathrm{}_{in}`$. The saturation magnetic moment $`m_{\mathrm{sat}}=I_\mathrm{s}V_\mathrm{d}/N_\mathrm{c}`$, where $`I_\mathrm{s}`$ is the saturation magnetization of the dot, plays role of the normalization factor in Eq.(2). Thus for the effective magnetic moments we have the bounding

$$m_{in,x}^2+m_{in,y}^21,$$
(3)
where $`x,`$ $`y`$ subscripts refer to the Cartesian components of the effective magnetic moments $`𝐦_{in}=m_{in,x}𝐞_x+m_{in,y}𝐞_y`$, where $`𝐞_x`$ and $`𝐞_y`$ are the Cartesian unit vectors. The ”softness” of $`𝐦_{in}`$ expressed by Eq.(3) is the important model aspect, which differs from the fundamental Brown’s postulate . The reason for this modification is that our model is formulated for sufficiently larger elements than classical micromagnetic approach. The total magnetization per dot per interaction center is given by

$$𝐦=\frac{1}{NN_\mathrm{c}}\underset{i=1}{\overset{N}{}}\underset{n=1}{\overset{N_\mathrm{c}}{}}𝐦_{in}.$$
(4)
The effective moments characterizing $`i`$th dot inside $`N`$-cluster are distributed around the dot center $`𝐑_i`$ and located at $`N_\mathrm{c}`$ positions

$$𝐗_{in}=𝐑_i+𝐫_n,n=1,2,\mathrm{},N_\mathrm{c},$$
(5)
where $`𝐫_n`$ are some relative coordinates of the interaction centers. From the assumption that identical dots are arranged into array it follows that system of $`𝐫_n`$ vectors is independent of the dot position inside the cluster.

In a quasi-two dimensional systems, where $`𝐦_{in}`$ is confined to $`xy`$ plane, reduced information about $`i`$th dot microstate is involved in $`2N_\mathrm{c}`$ dimensional row vector

$`\stackrel{~}{m}_i`$ $``$ $`[𝐦_{i1},𝐦_{i2},\mathrm{},𝐦_{iN_\mathrm{c}}]`$ (6)
$`=`$ $`[\underset{\mathrm{center}𝐗_{i1}}{\underset{}{m_{i1,x},m_{i1,y}}},\underset{\mathrm{center}𝐗_{i2}}{\underset{}{m_{i2,x},m_{i2,y}}},\mathrm{},\underset{\mathrm{center}𝐗_{iN_\mathrm{c}}}{\underset{}{m_{iN_\mathrm{c},x},m_{iN_\mathrm{c},y}}}].`$ (7)
We continue with the construction of relations associating the effective dot moments from Eq.(6) with the corresponding self-energies. Here the intradot self-energy is understood as a part of the total energy, which includes only anisotropy, exchange and intradot magnetostatic energy contributions. The Zeeman term and interdot magnetostatic terms, which do not contribute to the self-energy are defined independently of the neural network part of the model.

The construction of the self-energy formula $`E^{\mathrm{self}}(\stackrel{~}{m}_i)`$ is based on a proper choice of the set of special $`2N_\mathrm{c}`$-dimensional memorized vectors (the input patterns of neural networks) of the type Eq.(7). The memorized vectors are constructed by putting into row $`N_\mathrm{c}`$ two-dimensional vectors $`𝐩_n^{(q)}`$

$$\stackrel{~}{p}^{(q)}=[p_{1,x}^{(q)},p_{1,y}^{(q)},p_{2,x}^{(q)},p_{2,y}^{(q)},\mathrm{},p_{N_\mathrm{c},x}^{(q)},p_{N_\mathrm{c},y}^{(q)}].$$
(8)
Here, the superscript $`q`$ identifies so called feature . In our case it is an integer from the set $`\mathrm{\Lambda }_Q\{0,1,\mathrm{},Q1\}`$. In the analogy with Eq.(6), the subscripts of $`p_{n,x}^{(q)},p_{n,y}^{(q)}`$ run over the interaction centers $`𝐫_n`$.

To measure the differences between the configurations we have introduced the Euclidean norm

$$\stackrel{~}{m}=\stackrel{~}{m}\stackrel{~}{m}^\mathrm{T}$$
(9)
written here for some magnetic moment $`\stackrel{~}{m}`$ \[again encoded via the rule from Eq.(7)\]. The superscript $`\mathrm{T}`$ from Eq.(9) denotes the vector transposition.

We start the construction of $`E^{\mathrm{self}}(\stackrel{~}{m}_i)`$ by assuming that self-energy is known for $`Q`$ memorized vectors

$$E^{\mathrm{self}}(\stackrel{~}{p}^{(q)})=w^{(q)},q\mathrm{\Lambda }_Q,$$
(10)
where $`w^{(q)}`$ are free parameters of our model. The quality of the interpolation via RBFN depends on the choice of the basis functions and corresponding weights. Most convenient for our preliminary purposes seems to be the use of Nadaraya-Watson regression estimator . According to this, the self-energy input-output relation can be written

$$E^{\mathrm{self}}(\stackrel{~}{m})=\underset{q=0}{\overset{Q1}{}}w^{(q)}\psi ^{(q)}(\stackrel{~}{m}),$$
(11)
where

$$\psi ^{(q)}(\stackrel{~}{m})=\frac{\mathrm{exp}\left({\displaystyle \frac{Q}{d_{\mathrm{max}}^2}}\stackrel{~}{p}^{(q)}\stackrel{~}{m}^2\right)}{{\displaystyle \underset{q=0}{\overset{Q1}{}}}\mathrm{exp}\left({\displaystyle \frac{Q}{d_{\mathrm{max}}^2}}\stackrel{~}{p}^{(q)}\stackrel{~}{m}^2\right)}$$
(12)
are radial basis functions $`\{\psi ^{(q)}(\stackrel{~}{m})|`$, $`q\mathrm{\Lambda }_Q\}`$ satisfying the normalization conditions

$$\underset{q=0}{\overset{Q1}{}}\psi ^{(q)}(\stackrel{~}{m})=1,0\psi ^{(q)}(\stackrel{~}{m})1.$$
(13)
The dot index of $`\stackrel{~}{m}`$ was omitted whenever the distinctions between the individual dots is unimportant. The choice of the dispersion $`d_{\mathrm{max}}/(\sqrt{2}Q)`$ in Eq.(12), where

$`d_{\mathrm{max}}`$ $``$ $`\begin{array}{c}\\ \mathrm{max}\\ q,q^{}\mathrm{\Lambda }_Q\end{array}\stackrel{~}{p}^{(q)}\stackrel{~}{p}^{(q^{})}`$ (17)
is consistent with the recommendation . Speaking in terms of the neural networks $`\psi ^{(q)}(\stackrel{~}{m})`$ is the activation function, which determines $`q`$th neuron’s response to a given input $`\stackrel{~}{m}`$ and self-energy parameter $`w^{(q)},q\mathrm{\Lambda }_Q`$ is the optimized weight of the link between the input and output layer of the network.

Having established the formula for the calculation of the self-energy, it is easy to prove that if the systems of the memorized vectors is composed from the conjugate vector pairs $`\stackrel{~}{p}^{(q)}`$, $`\stackrel{~}{p}^{(q^{})}`$ with the same self-energy parameters

$$\stackrel{~}{p}^{(q)}=\stackrel{~}{p}^{(q^{})},w^{(q)}=w^{(q^{})}$$
(18)
RBFN self-energy form including the typical combination of terms

$$w^{(q)}\left(\mathrm{exp}((Q/d_{\mathrm{max}}^2)\stackrel{~}{p}^{(q)}\stackrel{~}{m}^2)+\mathrm{exp}((Q/d_{\mathrm{max}}^2)\stackrel{~}{p}^{(q)}+\stackrel{~}{m}^2)\right)$$

posses the reflection symmetry

$$E^{\mathrm{self}}(\stackrel{~}{m})=E^{\mathrm{self}}(\stackrel{~}{m}).$$
(19)
To analyze the configuration snapshots generated during the simulation process we have introduced so called feature map $``$. It associates any pattern vector $`\stackrel{~}{m}`$ with the feature $`q^{}\mathrm{\Lambda }_Q`$ . The feature $`q^{}`$ identifies the index of a nearest memorized vector $`\stackrel{~}{p}^{(q^{})}`$:

$`:q^{}`$ $`=`$ $`(\stackrel{~}{m}),`$ (20)
$`\stackrel{~}{p}^{(q^{})}\stackrel{~}{m}`$ $`=`$ $`\begin{array}{c}\\ \mathrm{min}\\ q\mathrm{\Lambda }_Q\end{array}\stackrel{~}{p}^{(q)}\stackrel{~}{m}.`$ (24)
The classification of the features performed for the whole cluster gives rise to the $`N`$component vector $`(\stackrel{~}{m}_1),(\stackrel{~}{m}_2),\mathrm{},(\stackrel{~}{m}_N)`$. The information about this vector can be concentrated to the form of the sample averages

$$\overline{n}(q)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\delta _{q,(\stackrel{~}{m}_i)},$$
(25)
where $`\delta `$ is the usual Kronecker symbol.

### 2.2 Interdot interactions and interactions with the external field

The standard assumption about the interdot interactions is that they are essentially magnetostatic . By using the concept of the effective moments and interaction centers one can construct easily the interdot interaction potential. For this aim we expressed the energy contribution $`E_{ik,j}^{\mathrm{c}\mathrm{d}}`$ consisting of the dipole-dipole interactions of the effective moment $`𝐦_{in}`$ \[located at the $`n`$th center of $`i`$th dot at the position $`𝐗_{in}`$\] with $`N_\mathrm{c}`$ effective moments of $`j`$th dot ($`ji`$)

$$E_{in,j}^{\mathrm{c}\mathrm{d}}=\lambda \underset{\alpha =x,y}{}m_{in,\alpha }\underset{\beta =x,y}{}\underset{s=1}{\overset{N_\mathrm{c}}{}}A_{injs,\alpha \beta }m_{js,\beta },$$
(26)
where

$`A_{injs,\alpha \beta }`$ $`=`$ $`J_{\alpha \beta }\left({\displaystyle \frac{𝐗_{in}𝐗_{js}}{a}}\right),J_{\alpha \beta }(𝐱)={\displaystyle \frac{|𝐱|^2\delta _{\alpha \beta }3x_\alpha x_\beta }{|𝐱|^5}}.`$ (27)
Because $`𝐗_{in}`$, $`𝐗_{js}`$ are scaled by the lattice spacing $`a`$, the energy dimension is absorbed into dipolar constant

$$\lambda =\frac{m_{\mathrm{sat}}^2}{4\pi \mu _0\mu _\mathrm{r}a^3},$$
(28)
where $`\mu _\mathrm{r}`$ is the relative permeability of the matrix. For the interdot magnetostatic energy of dot pair $`(i,j)`$ we obtain the expression

$$E_{ij}^{\mathrm{d}\mathrm{d}}=\underset{n=1}{\overset{N_\mathrm{c}}{}}E_{in,j}^{cd}.$$
(29)
The key remaining contribution is the Zeeman energy. For $`i`$th dot interacting with the external magnetic field $`𝐇`$ we obtain

$$E_i^\mathrm{H}=m_{\mathrm{sat}}𝐇\underset{n=1}{\overset{N_\mathrm{c}}{}}𝐦_{in}.$$
(30)
In further, to characterize the external field, we have used the reduced undimensional field $`𝐡=m_{\mathrm{sat}}𝐇/\lambda `$, its components $`h_x=𝐡𝐞_x`$ $`=h\mathrm{cos}\theta `$, $`h_y=𝐡𝐞_y`$$`=h\mathrm{sin}\theta `$ and polar angle $`\theta `$. The final form of total energy functional is then given by

$$E=\underset{i=1}{\overset{N}{}}\left\{E^{\mathrm{self}}(\stackrel{~}{m}_i)+E_i^\mathrm{H}+\underset{j=i+1}{\overset{N}{}}E_{ij}^{\mathrm{d}\mathrm{d}}\right\}.$$
(31)
The schematic view on its structure is displayed in Fig.1. According to the scheme of the computation, the interacting dots have their counterpart in the interacting RBFN blocks. This aspect makes our formulation close to the concept of the interacting neural networks .

## 3 The application to the array of rectangle ultrathin magnetic dots

In the concrete case we have modelled the ultrathin (quasi-two dimensional) rectangle magnetic dots of the square profile $`\mathrm{}_\mathrm{d}\times \mathrm{}_\mathrm{d}`$ and very small height $`(V_\mathrm{d}/\mathrm{}_\mathrm{d}^2)`$ $`\mathrm{}_\mathrm{d}`$. Each dot is replaced by the four moments $`N_\mathrm{c}=4`$ at the positions $`𝐫_n=(1/4)\mathrm{}_\mathrm{d}`$ $`[𝐞_x\mathrm{cos}\left(\frac{\pi }{2}(n1)\right)`$ $`+𝐞_y\mathrm{sin}\left(\frac{\pi }{2}(n1)\right)]`$ (see Fig.2(a)). The magnetic moment of $`\mathrm{}_{in}`$ element \[see Fig.3(b)\] maintains the saturation value $`m_{\mathrm{sat}}=I_\mathrm{s}V_\mathrm{d}/4`$. After this Eq.(28) modifies to the form $`\lambda =V_\mathrm{d}^2I_\mathrm{s}^2/(64\pi \mu _0\mu _\mathrm{r})`$.

We have proposed model including $`Q=11`$ memorized patterns (see Table 1) and four desired self-energy parameters

$`w^{(q)}=\{\begin{array}{ccc}_0\hfill & \text{for}\hfill & q=0,\hfill \\ _\mathrm{v}\hfill & \text{for}\hfill & q=1,2,\hfill \\ _\mathrm{p}\hfill & \text{for}\hfill & q=3,4,5,6,\hfill \\ _\mathrm{d}\hfill & \text{for}\hfill & q=7,8,9,10.\hfill \end{array}`$ (36)
The additional parameters of the model are reduction coefficients $`\kappa _\mathrm{v}1`$, $`\kappa _\mathrm{p}1`$, $`\kappa _\mathrm{d}1`$ (Table 1) introduced to modify the size of memorized effective moments. These coefficients describe the deviations of moments from the saturated value. According to Eq.(36) one can introduce four subsets of vectors $`\{\stackrel{~}{p}^{(0)}\}`$, $`\{\stackrel{~}{p}^{(1)},\stackrel{~}{p}^{(2)}\}`$, $`\{\stackrel{~}{p}^{(3)},\stackrel{~}{p}^{(4)},\stackrel{~}{p}^{(5)},\stackrel{~}{p}^{(6)}\}`$, $`\{\stackrel{~}{p}^{(7)},\stackrel{~}{p}^{(8)},\stackrel{~}{p}^{(9)},\stackrel{~}{p}^{(10)}\}`$. Within to each subset, the vectors correspond to the same self-energy. This system includes the vector pairs of the opposite sign $`p^{(2)}=p^{(1)}`$, $`p^{(5)}=p^{(3)}`$, $`p^{(6)}=p^{(4)}`$, $`p^{(9)}=p^{(7)}`$, $`p^{(10)}=p^{(8)}`$ and zero memorized vector $`\stackrel{~}{p}^{(0)}`$. This structure guarantee the reflection symmetry of the self-energy given by Eq.(19).

The exceptional vector $`\stackrel{~}{p}^{(0)}`$ concerns the integral information from multidomain or chaotic magnetization modes of the oscillatory or chaotic character . Its occurrence is a signature of uncertainty in description of a high momentum magnetization modes. Among the patterns memorized and restored by RBFN, we focussed attention to the vortex magnetization modes . The next two vectors $`\stackrel{~}{p}^{(1)}`$ ,$`\stackrel{~}{p}^{(2)}`$ encode the symmetric vortex and counter-vortex configurations in Fig.3(c), Table 1. Similarly, as in the case of $`\stackrel{~}{p}^{(0)}`$, the total magnetic moment of the symmetric memorized vortex is zero. For the rectangular ultrathin isolated dots of square profile and small crystalline anisotropy, the vortex type of magnetic ordering was revealed by the Monte-Carlo simulations . This finding was confirmed by the experiments . Vortex modes were also detected by the simulations on a cubic particles  for a weak or zero external magnetic fields. The system of vectors $`\{\stackrel{~}{p}^{(3)}`$, $`\stackrel{~}{p}^{(4)}`$, $`\stackrel{~}{p}^{(5)}`$, $`\stackrel{~}{p}^{(6)}\}`$ belonging to the Stoner-Wohlfart type of single-domain particle  is represented by the four parallel effective moments. This ordering can also occur by virtue of the external magnetic or magnetostatic fields. The remaining intradot configurations labeled by $`q=7,8,9,10`$ should be the potential sources of (shape) anisotropy. Let us to note that RBFN approach is not sensitive to the physical nature of the anisotropy.

The principal question arises how to determine seven single-dot parameters $`\kappa _\mathrm{v}`$, $`\kappa _\mathrm{p}`$, $`\kappa _\mathrm{d}`$, $`_0`$, $`_\mathrm{v}`$, $`_\mathrm{p}`$ and $`_\mathrm{d}`$. Further work is needed to combine the present simulations with the micromagnetic approaches (see e.g.) incorporating the algorithms of the neural network learning . In this paper the magnetic configurations were selected and parametrized in heuristic manner.

## 4 The implementation of simulated annealing method

Simulated annealing  is an optimization technique which operates in a manner analogous to the physical process of annealing. In this section we discuss some details of its implementation to dot array model. The subject of minimization is energy functional Eq.(31) of the effective magnetic moments. The main parts of adopted algorithm are:

* Initial state $`\stackrel{~}{m}_i(t=0)`$ is generated (or read from the data file).
* Cooling schedule. The pseudotemperature $`T(t)=T_0\mathrm{exp}(t/t_0))`$ relaxes as a function of the discrete time $`t=0,1,\mathrm{},t_{\mathrm{max}}`$.
* The update equation is based on the standard algorithm of Metropolis .
* Single moment moves. The lattice dot index $`i\{1,2,\mathrm{},N\}`$, and interaction center index $`n\{1,2,\mathrm{},N_\mathrm{c}\}`$ are chosen randomly. At given time $`t`$, the current effective moment $`𝐦_{in}(t)`$ undergoes the stochastic single-moment move $`𝐦_{in}^{\mathrm{trial}}(t+1)=𝐦_{in}(t)+v_\mathrm{m}𝐘(t),`$ where $`𝐘`$ is two-dimensional vector generated over uniform distribution in $`<1,1>\times <1,1>`$. The generation is repeated until the constrain $`|𝐦_{in}(t)+v_\mathrm{m}𝐘(t)|1`$ is satisfied; $`v_\mathrm{m}`$ is the parameter, which controls the remagnetization speed $`\left|𝐦(t+1)𝐦(t)\right|v_\mathrm{m}/(NN_\mathrm{c})`$. The encoding of $`𝐦_{in}^{\mathrm{trial}}(t)`$ gives rise to the row vector $`\stackrel{~}{m}_i^{\mathrm{trial}}`$(t). By taking into account Eq.(31) with the concrete form of the Nadaraya-Watson self-energy estimator \[Eq.(11)\] and expression for dipole-dipole interaction Eq.(26) we obtained for the energy variation $`\mathrm{\Delta }E_{in}`$ associated with the elementary move from $`𝐦_{in}(t)`$ to $`𝐦_{in}^{\mathrm{trial}}(t+1)`$ as

$`\mathrm{\Delta }E_{in}={\displaystyle \underset{q=0}{\overset{Q1}{}}}w^{(q)}\left[\psi ^{(q)}(\stackrel{~}{m}_i^{\mathrm{trial}})\psi ^{(q)}(\stackrel{~}{m}_i)\right]`$ (37)
$`\lambda {\displaystyle \underset{\alpha =x,y}{}}\left(m_{in,\alpha }^{\mathrm{trial}}m_{in,\alpha }\right)\left[h_\alpha {\displaystyle \underset{\beta =x,y}{}}{\displaystyle \underset{j=1,ji}{\overset{N}{}}}{\displaystyle \underset{s=1}{\overset{N_\mathrm{c}}{}}}A_{injs,\alpha \beta }m_{js,\beta }\right]`$
* Acceptance criteria If $`\mathrm{\Delta }E_{in}0`$, the trial configuration is accepted automatically and $`𝐦_{in}(t+1)=𝐦_{in}^{\mathrm{trial}}(t+1)`$. If $`\mathrm{\Delta }E_{in}>0`$, the trial configuration is accepted if the Boltzmann factor $`\mathrm{exp}(\mathrm{\Delta }E_{in}/T)`$ is larger or equal to the random number generated uniformly over $`<0,1>`$.
* Quasistatic simulations. In the zero temperature limit, the acceptance criteria reduces to the absolute acceptance if $`\mathrm{\Delta }E_{in}0`$. For this dynamics the energy is monotonically decreasing function of time. Subsequently, the stochastic motion through the phase space tends to the accessible basin of attraction. In the matastable state the motion gets stuck for a fixed external magnetic field. The sequence of the metastable configurations obtained for a gradually changing external magnetic field was used for the calculation of the quasistatic hysteresis loops.
* Complex intradot moves allows to enhance the effectiveness of annealing and left the metastable states. The update according to $`\mathbf{3.1}`$ is supplemented by :

+ Interchange moves changing the location of two intradot configurations: $`\stackrel{~}{m}_i^{\mathrm{trial}}=\stackrel{~}{m}_j`$, $`\stackrel{~}{m}_j^{\mathrm{trial}}=\stackrel{~}{m}_i`$ belonging to two randomly chosen dots $`i`$ and $`j`$.
    + Feature moves starting with the random choice of the dot index $`i`$ and its feature $`q\mathrm{\Lambda }_Q`$, $`q(\stackrel{~}{m}_i)`$. The suggested move is then given by $`\stackrel{~}{m}_i^{\mathrm{trial}}=\stackrel{~}{p}^{(q)}`$.
    + Reflection moves defined by $`\stackrel{~}{m}_i^{\mathrm{trial}}=\stackrel{~}{m}_i`$ are of the special importance since they conserve the self-energy \[see Eq.(19)\]. Due to this property the simulated system can overcome easily the self-energy barriers.

For the moves $`\mathbf{3.2}𝐚,𝐛,𝐜`$ the acceptance probability is given by

$$\mathrm{min}\{\mathrm{\hspace{0.17em}1},\mathrm{exp}\left(\underset{n=1}{\overset{N_\mathrm{c}}{}}\mathrm{\Delta }E_{in}\right)\}.$$
* Stopping criteria. For $`t<t_{\mathrm{max}}`$, the annealing process follows from the step $`\mathrm{𝟐}`$ with $`tt+1`$ for the move $`\mathbf{3.1}`$, or with $`tt+N_\mathrm{c}`$ in the case of the complex intradot moves $`\mathbf{3.2}𝐚,𝐛,𝐜`$.

## 5 Numerical simulations

For the model defined in Sect.3 we performed the numerical simulations. We studied finite disk-shaped cluster $`R^{\mathrm{cluster}}=8a`$ including $`N=193`$ dots of the size $`l_\mathrm{d}=a/8`$. We assume that each dot has $`N_\mathrm{c}=4`$ centers at the square lattice \[see Fig.2\]. The parameters of memorized magnetic configurations are $`\kappa _\mathrm{v}=\kappa _\mathrm{d}=0.95`$, $`\kappa _\mathrm{p}=1`$. Consequently, $`d_{\mathrm{max}}=4`$. The typical parameters of the simulated annealing have been $`v_\mathrm{m}=0.2`$, $`t_0=200NN_\mathrm{c}`$, $`t_{\mathrm{max}}=5t_0`$, and $`\lambda =T_010\lambda `$. For the quasistatic simulations we used $`t_{\mathrm{max}}=100NN_\mathrm{c}`$.

The initial simulations were performed for a zero external field and zero self-energy parameters. In this case the Monte-Carlo minimization of the energy leads to the rapid fall-off of the energy towards the non-colinear antiferromagnetic chains. The configuration is displayed in Fig.8(a). From this follows that central part of this cluster corresponds to the noncolinear antiferromagnetic phase in agreement with magnetic configuration obtained for a system of cylindrical dots  and truncated dipolar moments . At the same time the surface moments which tend to be parallel to the cluster surface exhibit some kind of the frustration . The annealing leads to the ground state estimate $`E=E^{\mathrm{AF}}=38.77\lambda N`$.

The previous value can be understood as the threshold for the competition between the interdot and intradot structures. The natural way of the stabilization of the intradot vortices is to make the parallel structures of the moments $`(q=3,4,5,6)`$ energetically unfavorable. In the next we will analyzed more restrictive choice: $`w^{(q1,2)}|E^{\mathrm{AF}}|`$ with the calibration condition $`w^{(q=1,2)}=_\mathrm{v}=0`$.

The previously simulated system was purely magnetostatic. We follow with the simulations of the opposite kind of systems, where interdot interactions have been completely neglected. For these systems we have constructed a quasi-static hysteresis loops. The results have been obtained for the several combinations of the self-energy parameters (comparable with $`|E^{\mathrm{AF}}|`$). They are presented in Figs.4(a)-(d). In Table 2 we list differences between the desired self-energies and outputs of Nadaraya-Watson estimator. Because our parameters are free, the inaccuracy stemming from Nadaraya-Watson formula has no principal significance for the quality of the result. In the situations, where precision of output becomes to be more relevant, the sofisticated RBFN learning is required . The interesting situation has occurred for the self-energy parameters

$`_0=20\lambda ,_\mathrm{v}=0,_\mathrm{p}=60\lambda ,_\mathrm{d}=40\lambda `$ (38)
corresponding to Fig.4(a). The field dependence of $`\overline{n}(q)`$ observed during the remagnetization process for $`q=1,2`$ confirms the vortex stabilization around $`h0`$ and zero remanence. Qualitativelly similar behavior with the vortex annihilation and formation was observed in the experimental study . The situation changes dramatically when the strong interdot interactions are taken into account. Fig.5 shows that their influence causes the non-zero ramanence due to suppression of vortices: $`\overline{n}(q)0.07`$ for $`q=1,2`$. More detailed analysis of this fact has revealed that surface vortices are more stable than vortices from the central zones of the cluster. Two next figures 6(a),(b) show how the anisotropy of the memorized configurations is reflected by the hysteresis loops of magnetization components $`m_x(h)`$, $`m_y(h)`$ constructed for $`\theta =20^o`$.

In the case with the non-zero self-energy parameters interact magnetostatically, the choice of the initial conditions of the simulated annealing becomes to be more complicated. Several final configurations obtained by the annealing process are displayed in Figs.8(b)-(i). The preliminary runs evolving from the initial random state get stuck in the local minimum $`E=6.4\lambda N`$. The configuration of this metastable state is displayed in Fig.8(b). The application of the feature map \[see Fig.8(c)\] reveals that the clustering of the dot states resembles the formation of the homogeneous domains in $`Q`$state Potts model . The lowest energy $`E=1.58\lambda N`$ was obtained for the system initialized from the vortex state. No essential differences between the pure vortex and mixed vortex-countervortex initial conditions were observed, contrary to our expectation evoked by the study of dipolar system . The question of the helicity will require more detailed investigation. Fig.8(d) shows that magnetostatic deformation of vortices is rather pronounced feature. After the deformation, four-moment vortices acquire the nonzero magnetic moments and resemble the fans or vortices with the non-central Néel-type core. Their analysis via $``$ map shows the mixing of vortex features $`q=1,2`$ \[ $`\overline{n}(1)`$ $``$ $`\overline{n}(2)`$ $`0.5`$\] and separation of the clusters with different vorticity \[see Fig.8(d). This aposteriory finding confirms that parameters of Eq.(38) are sufficient for the stabilization of the vortex ground state for $`h=0`$. Two configurations displayed in Figs.8(e),(f) were obtained for the external magnetic field. Their energies are listed in Table 3. For $`h=h_x=5`$ the field-deformed vortices resemble the fans ordered into the large-scale wave-like structures. The waves are better visible from Fig.8(h) showing the detailed snapshot, where intradot magnetic moments are averaged for each dot separatelly. The chaining of field-oriented phase is also visible. For $`h=h_x=10`$ Fig.8(i) demonstrates the formation of the clusters with $`q=6`$.

To characterize the anisotropy, we have studied the angular dependence of the magnetization for $`h=5`$ and varying $`\theta <0^o,180^o>`$. In the simulations we have distinguished between the clockwise and counter-clockwise field rotation directions. After the annealing starting from the purely vortex state $`\stackrel{~}{m}_{1iN}(t=0)=\stackrel{~}{p}^{(1)}`$, $`\theta =0^o`$ ($`\theta =180^o`$) we performed the series of quasistatic remagnetization steps at zero temperature. These simulations have revealed the hard axes $`𝐞_x`$, $`𝐞_y`$ and easy axes $`𝐞_x\pm 𝐞_y`$ \[see Fig.7\]. In addition, the model system exhibits the angular hysteresis. These results demonstrate how the outputs of RBFN mimic the biaxial anisotropy and how the anisotropy is reinforced by the magnetostatic couplings.

## 6 Conclusions

We believe that very general method we have introduced in this paper will be stimulating for the people working in the field of the micromagnetic simulations of the nanoscale systems. From the point of view of magnetostatics, direct model improvement is possible in many ways: a) near dot interactions can be taken into account more accuratelly by including the rectangle-rectangle magnetostatic terms; b) to speed-up the computations and to extend the system size one can use the hierarchical summation . We thing that more realistic simulations will be possible after the finding of a closer relationship between the neural networks and outputs of the standard micromagnetic approaches.

Acknowledgement: This work was supported by the grant no.1/6020/99 and by the Polish-Slovak international grant.

## 7 Figure Captions


text:
# Tracer diffusion in granular shear flows

## I Introduction

Granular systems under rapid flow conditions can be modeled as a fluid of inelastic hard spheres. In the simplest model the grains are taken to be smooth so that the inelasticity is characterized through a constant coefficient of normal restitution. The essential difference with respect to molecular fluids is the absence of energy conservation yielding modifications of the usual hydrodynamic equations. Due to the kinetic energy dissipation in collisions, energy must be externally injected to the granular gas in order to achieve a stationary state. In some experimental situations, the granular system is driven into flow by the presence of a shear field. In this case, a steady state is possible when the amount of energy supplied by shearing work is balanced by that lost due to the inelastic cooling. The study of the rheological properties of this steady shear flow state has received a great attention in the past years, especially in the case of one-component systems. However, much less is known in the more complicated case of multicomponent mixtures of grains.

An interesting problem is the analysis of diffusion in granular shear flows. The understanding of mass transport in granular systems is of practical interest since, for instance, powders must frequently be mixed togheter before any sort of processing can begin. The self-diffusion phenomenon in granular flows has been studied earlier. Experimental studies include both systems with macroscopic flows  and vertically vibrated systems. Complementary computer simulation studies have been also carried out, with special emphasis on the influence of the solid volume fraction on the diffusive motion of the grains. In general, all these previous studies were limited to observing the diffusion in only one direction, usually the direction parallel to the velocity gradient. However, due to the anisotropy induced in the system by the presence of shear flow, a diffusion tensor is required to describe the diffusion process instead of a single diffusion coefficient. To the best of my knowledge, the only attempt to measure the elements of this self-diffusion tensor has been made by Campbell. He measured these elements via molecular dynamics simulations by using both particle tracking and through velocity correlations. Both methods were found to agree with reasonable accuracy. Very recently, the self-diffusion coefficients have been experimentally measured in a granular system under Couette flow by employing image technology.

In the context of kinetic theory, studies on granular flows in mixtures are scarcer. Most of them  are based on a Navier-Stokes description and, therefore, they are restricted to small velocity gradients, which for the steady simple shear flow is equivalent to the low-dissipation limit. For this reason, the diffusion is only characterized by a single coefficient which is not affected by the presence of the shear field. In addition, although these studies permit in principle different temperatures for the two species, they assume equal partial granular temperatures $`T_i`$ in the quasielastic limit. Nevertheless, given the intrinsic connection between the shear rate and dissipation in this problem, energy nonequipartition is expected as the restitution coefficient decreases. As a matter of fact, some recent results obtained in molecular dynamics simulations of granular sheared mixtures as well as in real experiments of vibrated mixtures in three  and two dimensions  clearly show the breakdown of energy equipartition. This implies that the temperatures $`T_i`$ are different from the mixture temperature $`T`$. The consequences of this effect on the transport properties are in general significant, as has been recently found in the freely cooling case. In conclusion, a consistent theory describing diffusion in granular shear flows must take into account both the tensorial character of the mass transport as well as the possibility of temperature differences.

The aim of this paper is to get the diffusion tensor in a binary granular mixture under simple shear flow in the framework of the Boltzmann equation. Due to the complexity of the general problem, here we consider the special case in which one of the components (say for instance, the species 1) is present in tracer concentration. The tracer problem is more amenable to analytically treatment. First, the tracer particles are directly enslaved to the granular gas and, second there are fewer parameters. Therefore, in this situation one can assume that the velocity distribution function $`f_2`$ of the excess component (granular gas) obeys a (closed) nonlinear Boltzmann equation while the velocity distribution function $`f_1`$ of the tracer particles (impurities) satisfies a Boltzmann-Lorentz equation. The starting point of our study is a recent solution of the set of Boltzmann equations for a binary mixture of inelastic hard spheres under shear flow. The corresponding Boltzmann-Lorentz equation for the impurities is solved by means of a perturbative scheme in powers of the gradient of the mole fraction of tracer particles. The main feature of this expansion is that the reference state around which we perturb is not restricted to small values of the shear rate, which for the steady shear flow problem is equivalent to arbitrary degree of dissipation. In the first order of the expansion, the tracer diffusion tensor is identified from the mass flux. Explicit expressions for the nonzero elements of this tensor are obtained by using a first Sonine polynomial approximation. These elements are given in terms of the restitution coefficients and the parametes of the mixture (masses and sizes).

The plan of the paper is as follows. In Sec. II we introduce the set of coupled Boltzmann equations describing the mixture and state the problem we are interested in. The state of the mixture in the absence of diffusion is analyzed and in particular, the nonzero elements of the pressure tensor $`𝖯_2`$ of the granular gas are obtained in the leading Sonine approximation. Comparison of these results with previous theories is also presented. The Section ends studying the state of tracer particles with special emphasis in the evaluation of the specific dependence of the temperature ratio $`T_1/T_2`$ on restitution coefficients, mass ratio, and size ratio. Section III deals with the perturbation scheme used to solve the Boltzmann-Lorentz equation of the tracer particles when the diffusion takes place in the system. This Section contains the main results of the paper since we determine the tracer diffusion tensor in the first order of the expansion of the concentration gradient. The dependence of the nonzero elements of this tensor on the different parameters of the problem is illustrated in the three dimensional case, showing a good qualitative agreement with Campbell’s simulations . Finally, in Sec. IV we close the paper with some concluding remarks.

## II Description of the problem: Granular mixture in simple shear flow

We consider a granular binary mixture composed by smooth inelastic disks ($`d=2`$) or spheres ($`d=3`$) of masses $`m_1`$ and $`m_2`$ and diameters $`\sigma _1`$ and $`\sigma _2`$. Collisions between particles are inelastic and characterized by three constant restitution coefficients $`\alpha _{11}`$, $`\alpha _{22}`$, and $`\alpha _{12}=\alpha _{21}`$, where $`\alpha _{ij}1`$ refers to the restitution coefficient for collisions between particles of species $`i`$ and $`j`$. In the low-density regime, the one-particle velocity distribution functions $`f_i(𝐫,𝐯_1;t)`$ ($`i=1,2`$) obey the set of nonlinear Boltzmann kinetic equations:

$$\left(\frac{}{t}+𝐯_1\right)f_i(𝐫,𝐯_1;t)=\underset{j}{}J_{ij}\left[𝐯_1|f_i(t),f_j(t)\right],$$
(1)
where the Boltzmann collision operator $`J_{ij}\left[𝐯_1|f_i,f_j\right]`$ describing the scattering of pairs of particles is

$`J_{ij}\left[𝐯_1|f_i,f_j\right]`$ $`=`$ $`\sigma _{ij}^{d1}{\displaystyle 𝑑𝐯_2𝑑\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠_{12})(\widehat{𝝈}𝐠_{12})}`$ (3)
$`\times \left[\alpha _{ij}^2f_i(𝐫,𝐯_1^{},t)f_j(𝐫,𝐯_2^{},t)f_i(𝐫,𝐯_1,t)f_j(𝐫,𝐯_2,t)\right].`$
Here, $`d`$ is the dimensionality of the system, $`\sigma _{ij}=\left(\sigma _i+\sigma _j\right)/2`$, $`\widehat{𝝈}`$ is a unit vector along their line of centers, $`\mathrm{\Theta }`$ is the Heaviside step function and $`𝐠_{12}=𝐯_1𝐯_2`$. In addition, the primes on the velocities denote the initial values $`\{𝐯_1^{},𝐯_2^{}\}`$ that lead to $`\{𝐯_1,𝐯_2\}`$ following a binary collision:

$$𝐯_1^{}=𝐯_1\mu _{ji}\left(1+\alpha _{ij}^1\right)(\widehat{𝝈}𝐠_{12})\widehat{𝝈},𝐯_2^{}=𝐯_2+\mu _{ij}\left(1+\alpha _{ij}^1\right)(\widehat{𝝈}𝐠_{12})\widehat{𝝈},$$
(4)
where $`\mu _{ij}=m_i/\left(m_i+m_j\right)`$. At a hydrodynamic level, the relevant quantities are the number densities $`n_i`$, the flow velocity $`𝐮`$, and the “granular” temperature $`T`$. They are defined in terms of moments of the distribution $`f_i`$ as

$$n_i=𝑑𝐯f_i(𝐯),\rho 𝐮=\underset{i}{}\rho _i𝐮_i=\underset{i}{}𝑑𝐯m_i𝐯f_i(𝐯),$$
(5)
$$nT=\underset{i}{}n_iT_i=\underset{i}{}𝑑𝐯\frac{m_i}{d}V^2f_i(𝐯),$$
(6)
where $`n=n_1+n_2`$ is the total number density, $`\rho =\rho _1+\rho _2`$ is the total mass density, and $`𝐕=𝐯𝐮`$ is the peculiar velocity. Equations (5) and (6) also define the flow velocity $`𝐮_i`$ and the partial temperature $`T_i`$ of species $`i`$, which measures the mean kinetic energy of species $`i`$.

The collision operators conserve the particle number of each species and the total momentum, but the total energy is not conserved. This implies that

$$\underset{i,j}{}𝑑𝐯\frac{1}{2}m_iV^2J_{ij}[𝐯|f_i,f_j]=\frac{d}{2}nT\zeta ,$$
(7)
where $`\zeta `$ is identified as the “cooling rate” due to inelastic collisions among all species. At a kinetic level, it is also convenient to discuss energy transfer in terms of the “cooling rates” $`\zeta _i`$ for the partial temperatures $`T_i`$. They are defined as

$$\zeta _i=\frac{1}{dn_iT_i}\underset{j}{}𝑑𝐯m_iV^2J_{ij}[𝐯|f_i,f_j].$$
(8)
The total cooling rate $`\zeta `$ can be expressed in terms of $`\zeta _i`$ as

$$\zeta =T^1\underset{i}{}x_iT_i\zeta _i,$$
(9)
where $`x_i=n_i/n`$ is the mole fraction of species $`i`$.

Let us now describe the problem we are interested in. We consider a granular binary mixture in which the masses and sizes of both species are arbitrary. Our aim is to analyze a diffusion problem in a mixture of grains driven into flow by the action of shear forces. Here, we consider the simple case of the tracer limit, namely, a binary mixture in which one of the components (say, for instance, 1) is present in tracer concentration ($`x_11`$). In the tracer limit, one expects that the state of the granular gas is not affected by the presence of the tracer particles so that its velocity distribution function $`f_2`$ obeys a (closed) nonlinear Boltzmann equation. Furthermore, the mole fraction of tracer particles is so small that their mutual interactions can be neglected in the kinetic equation of $`f_1`$. As a consequence, the velocity distribution function of tracer particles $`f_1`$ satisfies a (linear) Boltzmann-Lorentz equation. Let us start describing the state of the mixture in the absence of diffusion.

### A Granular gas (excess component)

We assume that the granular gas is subjected to the simple shear flow. From a macroscopic point of view, this state is characterized by a constant linear velocity profile $`𝐮=𝐮_2=𝖺𝐫`$, where the elements of the tensor $`𝖺`$ are $`a_k\mathrm{}=a\delta _{kx}\delta _\mathrm{}y`$, $`a`$ being the constant shear rate. In addition, the partial density $`nn_2`$ and the granular temperature $`TT_2`$ are uniform. The temporal variation of the granular temperature $`T_2`$ can be obtained from the Boltzmann equation (1) as

$$\frac{p_2}{t}=aP_{2,xy}\frac{d}{2}\zeta _2p_2,$$
(10)
where $`p_2=n_2T_2`$,

$$𝖯_2=m_2𝑑\mathrm{𝐯𝐕𝐕}f_2$$
(11)
is the pressure tensor of the gas, and in the tracer limit

$$\zeta _2=\frac{1}{dn_2T_2}𝑑𝐯m_2V^2J_{22}[𝐯|f_2,f_2].$$
(12)
The balance of energy (10) shows the different nature of this state for molecular and granular systems. While for elastic fluids ($`\zeta _2=0`$) the temperature increases monotonically in time due to the viscous heating term $`aP_{2,xy}`$, a steady state is possible for granular systems when the viscous heating is exactly compensated by the collisional cooling term $`(d/2)p_2\zeta _2`$. In that case, the shear stress $`P_{2,xy}`$ and the cooling rate $`\zeta _2`$ are related by

$$aP_{2,xy}=\frac{d}{2}\zeta _2p_2.$$
(13)
As a consequence, for a given shear rate $`a`$, the (steady) temperature $`T_2`$ is a function of the restitution coefficient $`\alpha _{22}`$. This steady state is what we want to analyze.

The simple shear flow becomes spatially uniform when one refers the velocities of the particles to a frame moving with the flow velocity $`𝐮`$: $`f_2(𝐫,𝐯)f_2(𝐕)`$. Consequently, the stationary Boltzmann equation for the excess component becomes

$$aV_y\frac{}{V_x}f_2(𝐕)=J_{22}[𝐕|f_2,f_2].$$
(14)
We are mainly interested in computing the nonzero elements of the pressure tensor $`𝖯_2`$. These elements can be obtained by multiplying the Boltzmann equation (14) by $`m_2V_kV_{\mathrm{}}`$ and integrating over $`𝐕`$. The result is

$$a_{km}P_{2,\mathrm{}m}+a_\mathrm{}mP_{2,km}=A_k\mathrm{},$$
(15)
where

$$A_k\mathrm{}=m_2𝑑𝐕V_kV_{\mathrm{}}J_{22}[𝐕|f_2,f_2].$$
(16)
To get an explicit expression for $`𝖯_2`$ one needs to compute the right hand side of the set of equations (15). This requires the explicit knowledge of $`f_2`$, which is not known even in the elastic limit. However, one expects to get a good estimate of the low moments of the Boltzmann collisional operator by expanding $`f_2`$ in Sonine polynomials and then truncate the series after the first few terms. This approach is similar to the usual moment method for solving the Boltzmann equation in the elastic case. Therefore, we take the leading Sonine approximation:

$$f_2(𝐕)f_{2,M}(𝐕)\left[1+\frac{m_2}{2T_2}\left(\frac{P_{2,k\mathrm{}}}{p_2}\delta _k\mathrm{}\right)\left(V_kV_{\mathrm{}}\frac{1}{d}V^2\delta _k\mathrm{}\right)\right],$$
(17)
where $`f_{2,M}`$ is a Maxwellian distribution at the temperature of the gas, i.e.,

$$f_{2,M}(𝐕)=n_2\left(\frac{m_2}{2\pi T_2}\right)^{d/2}\mathrm{exp}\left(\frac{m_2V^2}{2T_2}\right).$$
(18)
With the approximation (17), the integrals appearing in the expressions of the cooling rate $`\zeta _2`$ and the collisional moment $`A_k\mathrm{}`$ can be explicitly evaluated. The details of the calculation are given in the Appendix. This allows us to get the explicit expressions of the non-zero elements of $`𝖯_2`$. They are given by

$$P_{2,yy}^{}=P_{2,zz}^{}=\mathrm{}=P_{2,dd}^{}=\frac{d+1+(d1)\alpha _{22}}{2d+33\alpha _{22}},$$
(19)
$$P_{2,xy}^{}=4d\frac{d+1+(d1)\alpha _{22}}{(1+\alpha _{22})(2d+33\alpha _{22})^2}a^{},$$
(20)
$$P_{2,xx}^{}=d(d1)P_{2,yy}^{},$$
(21)
$$a^2=\frac{d+2}{32d}\frac{(1+\alpha _{22})(2d+33\alpha _{22})^2(1\alpha _{22}^2)}{d+1+(d1)\alpha _{22}}.$$
(22)
Moreover, the (reduced) cooling rate $`\zeta _2^{}`$ is

$$\zeta _2^{}=\frac{d+2}{4d}(1\alpha _{22}^2).$$
(23)
Here, $`𝖯_2^{}=𝖯_2/p_2`$, $`a^{}=a/\nu `$, $`\zeta _2^{}=\zeta _2/\nu `$, and $`\nu `$ is a characteristic collision frequency given by $`\nu =p_2/\eta `$, with $`\eta `$ is the shear viscosity coefficient of the gas in the elastic limit, i.e.,

$$\eta =\frac{d+2}{8}\pi ^{(d1)/2}\mathrm{\Gamma }(d/2)\sigma _2^{(d1)}(m_2T_2)^{1/2}.$$
(24)
Equations (19)–(23) generalize previous results derived in the three dimensional case. The expression (22) clearly indicates the intrinsic connection between the velocity gradient and dissipation in the system. As a matter of fact, given that $`\alpha _{22}1`$, the range of (reduced) shear rates are defined in the interval $`0a^2(d+2)(3+2d)^2/32d(d+1)`$. The parameter $`a^{}`$ can be considered as the relevant nonequilibrium parameter of the system. In the elastic limit ($`\alpha _{22}=1`$ which implies $`a^{}=0`$), the equilibrium results of the molecular gas are recovered, i.e., $`P_{2,k\mathrm{}}^{}=\delta _k\mathrm{}`$. As said in the Introduction, the steady simple shear flow for a monocomponent granular fluid has been the subject of many previous works. Two interesting studies have been carried out by Jenkins and Richman for smooth inelastic disks  and by Brey et al. for a $`d`$ dimensional system. The latter description has been subsequently extended to dense gases . The approximated theory of Jenkins and Richman  is based in a generalized Maxwellian distribution to model the steady shear flow. On the other hand, Brey et al. solved a kinetic model equation of the Boltzmann equation and compared their predictions for the reduced temperature, shear stress, and normal stress differences with Monte Carlo simulations. Comparison between kinetic model results and simulation shows in general a good agreement. Our results differ from those reported in Refs.  and , although for practical purposes, the discrepancies between the different approaches are quite small, even for moderate values of $`\alpha _{22}`$. As an illustration, in Fig. 1 we compare the different theories for the (reduced) elements of the pressure tensor $`𝖯_2^{}`$ in a two-dimensional system ($`d=2`$). It is seen that the agreement is remarkable, although the discrepancies slightly increase as the restitution coefficient decreases.

### B Tracer component

In the absence of diffusion, the velocity distribution function $`f_1`$ of tracer particles satisfies the steady kinetic equation

$$aV_y\frac{}{V_x}f_1=J_{12}[𝐕|f_1,f_2].$$
(25)
The most relevant moment at this level of approximation is the pressure tensor $`𝖯_1`$ defined as

$$𝖯_1=m_1𝑑\mathrm{𝐕𝐕𝐕}f_1(𝐕),$$
(26)
The nonzero elements of $`𝖯_1`$ obey the set of equations:

$$a_{km}P_{1,\mathrm{}m}+a_\mathrm{}mP_{1,km}=B_k\mathrm{}$$
(27)
where

$$B_k\mathrm{}=m_1𝑑𝐕V_kV_{\mathrm{}}J_{12}[𝐕|f_1,f_2].$$
(28)
As done before in the case of $`𝖯_2`$, we estimate the collisional moment $`B_k\mathrm{}`$ by taking the leading Sonine approximation of $`f_1`$

$$f_1(𝐕)f_{1,M}(𝐕)\left[1+\frac{m_1}{2T_1}\left(\frac{P_{1,k\mathrm{}}}{p_1}\delta _k\mathrm{}\right)\left(V_kV_{\mathrm{}}\frac{1}{d}V^2\delta _k\mathrm{}\right)\right],$$
(29)
where $`p_1=n_1T_1`$ and now $`f_{1,M}`$ is a Maxwellian distribution at the temperature of the tracer particles $`T_1`$, i.e.,

$$f_{1,M}(𝐕)=n_1\left(\frac{m_1}{2\pi T_1}\right)^{d/2}\mathrm{exp}\left(\frac{m_1V^2}{2T_1}\right).$$
(30)
As will be shown later, the partial temperatures $`T_1`$ and $`T_2`$ are in general different so that the granular energy per particle is not equally distributed between both species.

Once the collisional moment $`B_k\mathrm{}`$ is determined (see the Appendix), the (reduced) non-zero components of $`𝖯_1^{}=𝖯_1/x_1p_2`$ can be easily obtained from Eq. (27). These components can be written in terms of the temperature ratio $`\gamma =T_1/T_2`$, the restitution coefficients $`\alpha _{22}`$ and $`\alpha _{12}`$, and the parameters of the mixture. After some algebra, one gets

$$P_{1,yy}^{}=P_{1,zz}^{}=\mathrm{}=P_{1,dd}^{}=\frac{F+HP_{2,yy}^{}}{G},$$
(31)
$$P_{1,xy}^{}=\frac{a^{}P_{1,yy}^{}HP_{2,xy}^{}}{G},$$
(32)
$$P_{1,xx}^{}=d\gamma (d1)P_{1,yy}^{},$$
(33)
where

$$F=\frac{\sqrt{2}}{2d}\left(\frac{\sigma _{12}}{\sigma _{22}}\right)^{d1}\mu _{12}\left(\frac{1+\theta }{\theta ^3}\right)^{1/2}(1+\alpha _{12})\left[1+\frac{\mu _{21}}{2}(d1)(1+\theta )(1+\alpha _{12})\right],$$
(34)
$`G`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{4d}}\left({\displaystyle \frac{\sigma _{12}}{\sigma _{22}}}\right)^{d1}\mu _{21}\left({\displaystyle \frac{1}{\theta (1+\theta )}}\right)^{1/2}(1+\alpha _{12})`$ (36)
$`\times \left\{2[(d+2)\theta +d+3]3\mu _{21}(1+\theta )(1+\alpha _{12})\right\},`$
$$H=\frac{\sqrt{2}}{4d}\left(\frac{\sigma _{12}}{\sigma _{22}}\right)^{d1}\mu _{12}\left(\frac{1}{\theta (1+\theta )}\right)^{1/2}(1+\alpha _{12})\left[3\mu _{21}(1+\theta )(1+\alpha _{12})2\right].$$
(37)
Here, $`\theta =m_1T_2/m_2T_1`$ is the mean square velocity of the gas particles relative to that of the tracer particles. To close the problem at this stage of approximation, it still remains to get the temperature ratio $`\gamma `$. It can be obtained for instance, from the requirements (13) and its corresponding counterpart for species 1. This yields

$$\gamma =\frac{\zeta _2^{}P_{1,xy}^{}}{\zeta _1^{}P_{2,xy}^{}},$$
(38)
where the cooling rate $`\zeta _1^{}=\zeta _1/\nu `$ for the tracer particles is (see the Appendix)

$$\zeta _1^{}=\frac{(d+2)\sqrt{2}}{4d}\left(\frac{\sigma _{12}}{\sigma _2}\right)^{d1}\mu _{21}\left(\frac{1+\theta }{\theta }\right)^{1/2}(1+\alpha _{12})\left[2\mu _{21}(1+\theta )(1+\alpha _{12})\right].$$
(39)
The solution to Eq. (38) gives $`\gamma `$ as a function of the restitution coefficients $`\alpha _{22}`$ and $`\alpha _{12}`$ and the mechanical parameters of the mixture, i.e., the mass ratio $`\mu =m_1/m_2`$ and the size ratio $`w=\sigma _1/\sigma _2`$. Except for some limiting cases, Eq. (38) must be solved numerically. Thus, in the elastic case ($`\alpha _{22}=\alpha _{12}=1`$), we recover the well-known equilibrium results with $`\gamma =1`$ and $`\theta m_1/m_2`$ as required by the equipartition theorem. In the case of mechanically equivalent particles ($`m_1=m_2`$, $`\alpha _{22}=\alpha _{12}`$, $`\sigma _1=\sigma _2`$), Eqs. (19)–(23) and (31)–(37) lead to $`𝖯_2^{}=𝖯_1^{}`$ and $`\zeta _2^{}=\zeta _1^{}`$, so that $`\gamma =1`$. Beyond the above limit cases, as expected, our results yield $`\gamma 1`$. The violation of energy equipartition in driven granular mixtures has been recently observed in molecular dynamics simulations of sheared mixtures  and in real experiments. This effect is generic of multicomponent granular systems and is consistent with previous results derived in the unforced case. To the best of my knowledge, the only previous theories including temperature differences have been proposed by Jenkins and Mancini  and by Huilin et al. . However, both works are phenomenological with no attempt to solve the kinetic equation. Instead, they assume that the velocity distribution function is a local Maxwellian. This is reasonable for estimating the dense gas collisional transfer contributions to the fluxes, but not for evaluating their kinetic contributions. Both theories are applicable to a general flow field. In particular, the results obtained by Jenkins and Mancini  for the temperature ratio in the low density limit for inelastic disks can be written as

$$\gamma =1\frac{1+\mu }{1+w}\left[\frac{1+w}{2\mu }(1\alpha _{12})\sqrt{\frac{\mu (1+\mu )}{2}}(1\alpha _{22})\right].$$
(40)
In Fig. 2 we plot the temperature ratio $`\gamma `$ versus the restitution coefficient $`\alpha `$ for a size ratio $`w=2`$ and three values of the mass ratio $`\mu `$ in the two-dimensional case ($`d=2`$). For the sake of simplicity, henceforth we will assume that the spheres or disks are made of the same material, i.e., $`\alpha \alpha _{22}=\alpha _{12}`$. Also for comparison, we show the prediction given by Eq. (40) in the case $`\mu =5`$. Important discrepancies between both theories appear even for values of $`\alpha `$ close to $`1`$. As a matter of fact, the theory of Jenkins and Mancini predicts a violation of energy equipartition much more significant than our theory. It must be remarked that the quantitative predictions of our theory at the level of the temperature ratio have been recently confirmed by Monte Carlo simulations. Regarding the influence of the parameters of the mixture, we observe that for large mass ratios the temperature differences are quite significant, even for moderate dissipation (say $`\alpha 0.9`$). The temperature of the tracer particles is larger than that of the excess species when the tracer grains are heavier than the grains of the gas. This behavior has been also found in the recent computer simulations carried out in rapid shear flow .

## III Tracer diffusion under simple shear flow

We want to study the diffusion of tracer particles immersed in a bath (granular gas) subjected to the simple shear flow. The diffusion process is induced in the system by a weak concentration gradient $`x_1`$. However, given that the strength of the shear rate $`a`$ is arbitrary, the mass flux (which is generated by the gradient $`x_1`$) can be modified by the presence of the shear flow. As stated above, in the tracer limit the state of particles of species 1 is mainly governed by the collisions with particles of species 2, so that the self-collisions among particles 1 can be neglected. Thus, the kinetic equation governing the evolution of the velocity distribution function $`f_1`$ reads

$$\frac{}{t}f_1aV_y\frac{}{V_x}f_1+(V_k+a_k\mathrm{}r_{\mathrm{}})\frac{}{r_k}f_1=J_{12}[𝐕|f_1,f_2],$$
(41)
where here the derivative $`/r_k`$ is taken at constant $`𝐕`$. Tracer particles may freely exchange momentum and energy with the particles of the granular gas and, therefore, these are not invariants of the collision operator $`J_{12}[f_1,f_2]`$. Only the number density of tracer particles is conserved. More specifically, the mole fraction $`x_1`$ obeys the conservation law

$$\left(\frac{}{t}+a_k\mathrm{}r_{\mathrm{}}\frac{}{r_k}\right)x_1+\frac{𝐣_1}{m_1}=0,$$
(42)
where the mass flux $`𝐣_1`$ is defined as

$$𝐣_1=m_1𝑑\mathrm{𝐕𝐕}f_1(𝐕).$$
(43)
When $`\alpha _{22}=\alpha _{12}=1`$ (which gives $`a^{}=0`$), the well-known Fick law establishes a linear relationship between the mass flux $`𝐣_1`$ and the concentration gradient $`x_1`$. This law defines the diffusion coefficient. For finite values of the (reduced) shear rate $`a/\nu `$ (which means $`\alpha _{22}1`$), one expects that a generalized Ficks’s law holds but now a diffusion tensor rather than a scalar should appear. Our aim is to get this tensor in terms of $`\alpha _{22}`$, $`\alpha _{12}`$, $`\mu `$, and $`w`$. To this end, and assuming that the mole fraction $`x_1`$ is slightly nonuniform, we solve Eq. (41) by means of a perturbation expansion around a nonequilibrium state with arbitrary shear rate, which is equivalent to strong dissipation in the simple shear flow \[see Eq. (22)\]. Thus, we write

$$f_1=f_1^{(0)}+f_1^{(1)}+\mathrm{},$$
(44)
where $`f_1^{(k)}`$ is of order $`k`$ in $`x_1`$ but applies for arbitrary degree of dissipation since this distribution retains all the orders in $`a`$. The solution (44) qualifies as a normal solution since all the space and time dependence of $`f_1`$ occurs entirely through $`x_1(𝐫;t)`$ and their gradients. The zeroth-order approximation $`f_1^{(0)}`$ corresponds to the simple shear flow distribution but taking into account now the local dependence on the mole fraction $`x_1`$. Although the explicit form of $`f_1^{(0)}`$ is not exactly known, only the knowledge of its second-degree moments (related to the pressure tensor $`𝖯_1`$) is necessary to get the diffusion tensor in the Sonine approximation.

The kinetic equation for $`f_1^{(1)}`$ can be obtained from the Boltzmann-Lorentz equation (41) by collecting all the terms of first order in $`x_1`$:

$$\frac{}{t}f_1^{(0)}aV_y\frac{}{V_x}f_1^{(1)}+(V_i+a_k\mathrm{}r_{\mathrm{}})\frac{}{r_k}f_1^{(0)}=J_{12}[𝐕|f_1^{(1)},f_2].$$
(45)
According to the balance equation (42), one has

$$\frac{f_1^{(0)}}{t}=\frac{f_1^{(0)}}{x_1}\frac{x_1}{t}=a_k\mathrm{}r_{\mathrm{}}\frac{x_1}{r_k}\frac{f_1^{(0)}}{x_1}$$
(46)
where use has been made of the fact that the zeroth-order approximation to the mass flux vanishes, i.e., $`𝐣_1^{(0)}=\mathrm{𝟎}`$. Moreover,

$$\frac{f_1^{(0)}}{r_k}=\frac{f_1^{(0)}}{x_1}\frac{x_1}{r_k}$$
(47)
Using (46) and (47), Eq. (45) can be written as

$$\left(aV_y\frac{}{V_x}+\mathrm{\Lambda }\right)f_1^{(1)}=\frac{f_1^{(0)}}{x_1}\left(𝐕x_1\right),$$
(48)
where $`\mathrm{\Lambda }`$ is the Boltzmann-Lorentz collision operator

$$\mathrm{\Lambda }f_1^{(1)}=J_{12}[𝐕|f_1^{(1)},f_2].$$
(49)
The solution to Eq. (48) is proportional to $`x_1`$, i.e., it has the form

$$f_1^{(1)}(𝐕)=𝒜(𝐕)x_1$$
(50)
Substitution of this into Eq. (48) yields

$$\left(aV_y\frac{}{V_x}+\mathrm{\Lambda }\right)𝒜=\frac{f_1^{(0)}}{x_1}𝐕.$$
(51)
The first order approximation to the mass flux is given by

$`𝐣_1^{(1)}`$ $`=`$ $`m_1{\displaystyle 𝑑\mathrm{𝐕𝐕}f_1^{(1)}(𝐕)}`$ (52)
$`=`$ $`𝖣x_1,`$ (53)
where the tracer diffusion tensor is

$$D_k\mathrm{}=m_1𝑑𝐕V_k𝒜_{\mathrm{}}(𝐕).$$
(54)
The solution to the integral equation (51) allows one to determine the quantity $`𝒜`$. From this solution one can determine the tracer diffusion tensor by means of Eq. (54).

In order to get an explicit expression for the tensor $`𝖣`$ we need to know the quantity $`𝒜`$. A good estimate of $`𝒜`$ to evaluate the mass flux $`𝐣_1^{(1)}`$ is given by the first Sonine approximation, in which only the leading term in the expansion of $`𝒜(𝐕)`$ in Sonine polynomials is kept. Thus, we take the following approximation to $`𝒜`$:

$$𝒜(𝐕)\frac{1}{n_1T_1}𝐕𝖣f_{1,M}(𝐕).$$
(55)
Using Eq. (55), an equation for the diffusion tensor is easily derived from Eq. (51). The result can be written as

$$\left(𝖺+𝛀\right)𝖣=p_2𝖯_1^{},$$
(56)
where the nonzero components of $`𝖯_1^{}`$ are given by Eqs. (31)–(33) and we have introduced the tensorial quantity

$$𝛀=\frac{m_1}{n_1T_1}𝑑𝐕𝐕\mathrm{\Lambda }𝐕f_{1,M}.$$
(57)
The expresion of the tensor $`𝛀`$ has been obtained in the Appendix with the result

$$𝛀=\frac{2}{d}\frac{\pi ^{(d1)/2}}{\mathrm{\Gamma }(d/2)}n_2\mu _{21}\sigma _{12}^{d1}(2T_2/m_2)^{1/2}(1+\alpha _{12})\left[(1+\theta )\theta \right]^{1/2}\left[\left(1+\frac{d+1}{d+2}\theta \right)𝟙+\frac{\theta }{𝕕+\mathrm{𝟚}}𝖯_\mathrm{𝟚}^{}\right].$$
(58)
The solution of Eq. (56) is

$$𝖣=p_2\left(𝖺+𝛀\right)^1𝖯_1^{}$$
(59)
Equation (59) is the primary result of this paper. It provides an explicit expression of the tracer diffusion tensor of a granular binary mixture in simple shear flow. The elements of this tensor give all the information on the physical mechanisms involved in the diffusion of tracer particles in a strongly sheared granular gas. In the absence of shear rate (what is equivalent to $`\alpha _{22}=\alpha _{12}=1`$), $`D_k\mathrm{}=D_0\delta _k\mathrm{}`$, where

$$D_0=\frac{d}{4\sqrt{2}}\frac{\mathrm{\Gamma }(d/2)}{\pi ^{(d1)/2}\sigma _{12}^{d1}}\sqrt{\mu (1+\mu )}\left(m_2T_2\right)^{1/2}$$
(60)
is the tracer diffusion coefficient of a molecular gas. As the restitution coefficients decrease, rheological effects become important and the elements of the diffusion tensor are different from the one obtained in the equilibrium case. The dependence of the diffusion coefficients on the restitution coefficients $`\alpha _{22}`$ and $`\alpha _{12}`$ as well as on the mass ratio $`\mu `$ and the size ratio $`w`$ is highly nonlinear. As happens for elastic fluids, Eq. (59) shows that diffusion under simple shear flow is a very complex problem due basically to the anisotropy induced in the system by the shear flow.

To illustrate the dependence of the elements $`D_k\mathrm{}`$ on the parameters of the problem, let us consider a mixture of inelastic hard-spheres ($`d=3`$). According to Eq. (59), $`D_{xz}=D_{zx}=D_{yz}=D_{zy}=0`$, in agreement with the symmetry of the problem. Consequently, there are five relevant elements: the three diagonal and two ($`D_{xy}`$, $`D_{yx}`$) off-diagonal elements. In general, $`D_{xx}D_{yy}D_{zz}`$ and $`D_{xy}D_{yx}`$. The off-diagonal elements measure cross effects in the diffusion of particles induced by the shear flow. Thus, for instance $`D_{xy}`$ gives the transport of mass along the direction of the flow of the system ($`x`$ axis) due to a concentration gradient parallel to the gradient of the flow velocity ($`y`$ axis). Both off-diagonal elements are negative.

Before analyzing the influence of the mechanical parameters of the mixture on diffusion, it is instructive to explore the particular case of self-diffusion, i.e., when the tracer particles are mechanically equivalent to the gas particles. This situation involves only single-particle motion and it is therefore somewhat simpler to compute the diffusion coefficients. In particular, the temperature of the tracer particles is the same as that of the gas particles and so $`\theta =1`$ in Eqs. (31)–(37) and (58). In Fig. 3, we plot $`D_{xx}^{}D_{yy}^{}`$, $`D_{zz}^{}D_{yy}^{}`$, $`(D_{xx}^{}+D_{yy}^{}+D_{zz}^{})/3(1/3)D_{kk}^{}`$, $`D_{xy}^{}`$, and $`D_{yx}^{}`$ as functions of the restitution coefficient $`\alpha \alpha _{12}=\alpha _{22}`$. Here, $`D_{ij}^{}D_{ij}/D_0`$, with $`D_0`$ given by Eq. (60). We see that the deviation from the functional form for elastic collisions is quite important even for moderate dissipation. Thus, for instance at $`\alpha =0.8`$, $`D_{xx}^{}D_{yy}^{}0.76`$, $`D_{zz}^{}D_{yy}^{}0.046`$, $`(1/3)D_{kk}^{}1.18`$, $`D_{xy}^{}1.039`$, and $`D_{yx}^{}0.42`$. The figure also shows that the anisotropy of the system, as measured by the differences $`D_{xx}^{}D_{yy}^{}`$ and $`D_{zz}^{}D_{yy}^{}`$, grows with the inelasticity. This anisotropy is much more important in the plane of shear flow ($`D_{xx}^{}D_{yy}^{}`$) that in the plane perpendicular to the flow velocity ($`D_{zz}^{}D_{yy}^{}`$). This is basically due to the fact that $`P_{s,xx}^{}P_{s,yy}^{}=P_{s,zz}^{}`$ with $`s=1,2`$.

As said before, Campbell  has carried out molecular dynamics simulations to measure the nonzero elements of the self-diffusion tensor. In his work, the self-diffusion coefficients were nondimensionalized by the product of the shear rate and the particle diameter. In our units, this corresponds to the reduced tensor $`\stackrel{~}{D}_{ij}=D_{ij}^{}/a^{}`$. Although the solid fractions analyzed in his simulations prevent us in general from making a quantitative comparison between our theory (restricted to dilute gases) and his computer simulations, we observe that the general qualitative dependence of the self-diffusion tensor on dissipation agrees well with our results, at least for the lowest solid fraction considered. Thus, theory and simulation predict that the magnitude of the normal diffusion coefficients follow the pattern $`\stackrel{~}{D}_{xx}>\stackrel{~}{D}_{zz}>\stackrel{~}{D}_{yy}`$ while, in general, the elements $`\stackrel{~}{D}_{ij}`$ decrease as the restitution coefficient decrease. An exception to the latter rule is the element $`\stackrel{~}{D}_{xx}`$, which does not depend sensitively on $`\alpha `$. On the other hand, in Campbell’s simulation work , he found that the values of $`\stackrel{~}{D}_{xy}`$ were roughly of the same magnitude as $`\stackrel{~}{D}_{yy}`$ provided that the solid fraction is smaller than 0.4. This trend is not completely followed by our theory since the values of $`\stackrel{~}{D}_{xy}`$ and $`\stackrel{~}{D}_{yy}`$ are significantly different for highly inelastic spheres. Thus, for instance, at $`\alpha =0.8`$, $`\stackrel{~}{D}_{xy}2.43`$ and $`\stackrel{~}{D}_{yy}2.15`$ but $`\stackrel{~}{D}_{xy}2.54`$ and $`\stackrel{~}{D}_{yy}1.10`$ at $`\alpha =0.4`$.

The dependence of the diffusion coefficients $`D_{ij}^{}`$ on the restitution coefficient for different values of the mass ratio $`\mu `$ is illustrated in Fig. 4. In this case,we take a size ratio $`w=2`$ and two values of the mass ratio: $`\mu =2`$ and $`\mu =4`$. For a given value of the inelasticity, we observe that the deviations from the elastic results are more important as the tracer particles are heavier than the gas particles.

## IV Concluding remarks

In this paper, we have described diffusion of tracer particles in a granular gas subjected to the simple shear flow. We have been interested in the steady state where the effect of viscosity is compensated for by the dissipation in collisions. Under these conditions, the resulting diffusion is anisotropic and, thus, cannot be described by a single diffusion coefficient. Instead, it must be described by a diffusion tensor whose explicit determination has been the main objective of this work. In order to capture the essential aspects of such a nonlinear problem, we have considered here a granular mixture in the low density regime as a prototype granular system, which lends itself to a detailed description by means of the nonlinear Boltzmann kinetic equation.

We have been concerned with the physical situation where a weak concentration gradient coexists with a strong shear rate, which for the steady simple shear flow problem means strong dissipation. For this reason, the Boltzmann-Lorentz equation corresponding to the tracer particles has been solved by means of a perturbation expansion around a nonequilibrium sheared state. This implies that the different approximations of this expansion are nonlinear functions of the restitution coefficients as well as of the parameters of the mixture (mass and size ratios). To get explicit results, we have used a first Sonine polynomial approximation to evaluate the cooling rates and the collisional moments of the Boltzmann operators. The reliablity of this approximation has been recently assessed in the (pure) shear flow problem where it has been shown to agree very well with Monte Carlo simulations in the case of hard spheres .

The kinetic theory results show that the elements $`D_{ij}`$ of the diffusion tensor present a complex dependence on the restitution coefficients $`\alpha _{22}`$ and $`\alpha _{12}`$ and on the mass ratio $`\mu =m_1/m_2`$ and the size ratio $`w=\sigma _1/\sigma _2`$. In the elastic case, $`\alpha _{ij}=1`$, $`D_{ij}=D_0\delta _{ij}`$, where $`D_0`$ is given by Eq. (60) and one recovers the expression of the diffusion coefficient for normal fluids. The deviations of the tensor $`D_{ij}`$ from the scalar $`D_0`$ have two distinct origins. First, the presence of shear flow gives rise to the new tensorial term $`𝖺+𝛀`$ on the left hand side of Eq. (56) instead of the corresponding collision frequency of the elastic diffusion problem. Second, given that the tracer and fluid particles are mechanically different, the reference state (zeroth-order approximation of the expansion) of tracer particles is completely different from that of the gas particles. In particular, when $`\mu 1`$ and/or $`w1`$, the temperature ratio $`\gamma `$ is clearly different from $`1`$ (as can be seen in Fig. 2), confirming the breakdown of the energy equipartition. The effect of different temperatures for the tracer and gas particles is expressed by the appearance of $`\theta =\mu /\gamma `$ in $`𝖯_1`$ \[cf. Eqs. (31)–(33)\] and in $`𝛀`$. Each one of the two afforementioned effects is a different reflection of dissipation present in the system.

A simple case is the self-diffusion problem, i.e., when the tracer and gas particles are mechanically equivalent. In this case, $`𝖯_2^{}=𝖯_1^{}`$ and $`\gamma =1`$. This situation has been previously studied by Campbell  by means of molecular dynamics simulations. As has been discussed in Sec. III, our predictions for the self-diffusion tensor agree qualitatively well with these simulations. On the other hand, when the tracer and gas particles are mechanically different, to my knowledge no previous studies on the diffusion tensor under shear flow have been made. As pointed out in the Introduction, most of the works on granular mixtures  are based on the Chapman-Enskog expansion around a local equilibrium state up to the Navier-Stokes order, and therefore, they are restricted to the low-dissipation limit in the simple shear flow. In addition, they also assume a single temperature to describe the mixture. A more careful calculation which takes into account temperature differences has been recently made by Garzó and Dufty. They have obtained explicit expressions for the transport coefficients of a granular binary mixture in terms of the restitution coefficients and the parameters of the mixture. Since these results have been derived taking the freely cooling state as the reference one, the diffusion is characterized by a single scalar coefficient that cannot be directly compared to the diffusion tensor obtained here. However, it would be interesting to compare the diffusion results obtained here in the driven sheared case with those found in the unforced case. In Fig. 5, we have compared the behavior of the scalar $`\frac{1}{3}D_{kk}^{}`$ (which can be understood as a generalized mutual diffusion coefficient in a sheared mixture) with the (reduced) diffusion coefficient $`D^{}`$ obtained in Ref.  in the tracer limit ($`x_10`$). We observe that, although the reference states in both descriptions are very different, the dependence of both diffusion coefficients on dissipation is quite similar since they increase as $`\alpha `$ increases. This trend is more significant in the unforced case than in the sheared case.

The evaluation of the diffusion tensor for practical purposes requires the truncation of a Sonine polynomial expansion. In the case of the (pure) simple shear flow problem, recent Monte Carlo simulations  have shown the accuracy of the leading order truncation. We expect that this agreement may be extended to the elements of the diffusion tensor for a wide range of values of dissipation. Exceptions to this agreement could be disparate mass binary mixtures (e.g., electron-proton systems) for which the first Sonine solution is not perhaps a good approximation and higher-order terms should be considered. We hope that the results derived in this paper will stimulate the performance of computer simulations to check the quality of the approximations given here for the diffusion tensor. Given the difficulties associated with molecular dynamics simulations in the low-density regime, one could perhaps use the direct simulation Monte Carlo method  which is being shown to be fruitful in the context of granular systems.

###### Acknowledgements.

Partial support from the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM2001-0718 is ackowledged. The author is grateful to Andrés Santos for suggestions for the improvement of the manuscript.

## A Evaluation of $`\zeta _i`$, $`𝖠`$, $`𝖡`$, and $`𝛀`$

In this Appendix we evaluate the cooling rates $`\zeta _i`$, the collisional velocity moments $`𝖠`$, and $`𝖡`$, and the tensorial collision frequency $`𝛀`$ by using the corresponding leading Sonine approximations.

### 1 Evaluation of $`\zeta _i`$

The cooling rate $`\zeta _1`$ of the tracer particles is defined as

$$\zeta _1=\frac{1}{dn_1T_1}𝑑𝐕_1m_1V_1^2J_{12}[𝐕_1|f_1,f_2].$$
(A1)
A useful identity for an arbitrary function $`h(𝐕_1)`$ is given by

$`{\displaystyle 𝑑𝐕_1h(𝐕_1)J_{12}\left[𝐕_1|f_1^{(0)},f_2\right]}`$ $`=`$ $`\sigma _{12}^{d1}{\displaystyle 𝑑𝐕_1𝑑𝐕_2f_1^{(0)}(𝐕_1)f_2(𝐕_2)}`$ (A3)
$`\times {\displaystyle }d\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠)(\widehat{𝝈}𝐠)[h(𝐕_1^{\prime \prime })h(𝐕_1)],`$
with

$$𝐕_1^{^{\prime \prime }}=𝐕_1\mu _{21}(1+\alpha _{12})(\widehat{𝝈}𝐠)\widehat{𝝈}.$$
(A4)
Using (A3), Eq. (A1) can be written as

$`\zeta _1`$ $`=`$ $`{\displaystyle \frac{m_1}{dn_1T_1}}\sigma _{12}^{d1}(1+\alpha _{12})\mu _{21}{\displaystyle 𝑑𝐕_1𝑑𝐕_2f_1(𝐕_1)f_2(𝐕_2)𝑑\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠_{12})(\widehat{𝝈}𝐠_{12})^2}`$ (A6)
$`\times \left[\mu _{21}(1+\alpha _{12})(\widehat{𝝈}𝐠_{12})2(𝐕_1\widehat{𝝈})\right]`$
To perform the angular integrations, we need the results

$$𝑑\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠_{12})(\widehat{𝝈}𝐠_{12})^n=\beta _ng_{12}^n,$$
(A7)
$$𝑑\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠_{12})(\widehat{𝝈}𝐠_{12})^n\widehat{𝝈}=\beta _{n+1}g_{12}^{n1}𝐠_{12},$$
(A8)
where

$$\beta _n=\pi ^{(d1)/2}\frac{\mathrm{\Gamma }\left((n+1)/2\right)}{\mathrm{\Gamma }\left((n+d)/2\right)}.$$
(A9)
Thus, the integration over $`\widehat{𝝈}`$ in Eq. (A6) leads to

$`\zeta _1`$ $`=`$ $`{\displaystyle \frac{m_1}{dn_1T_1}}\sigma _{12}^{d1}\beta _3(1+\alpha _{12})\mu _{21}{\displaystyle 𝑑𝐕_1𝑑𝐕_2f_1(𝐕_1)f_2(𝐕_2)}`$ (A11)
$`\times \left[\mu _{21}g_{12}^3(1\alpha _{12})+2g_{12}(𝐠_{12}𝐆_{12})\right],`$
where $`𝐆_{12}=\mu _{12}𝐕_1+\mu _{21}𝐕_2`$ is the center of mass velocity. Now, we take the Sonine approximations given by (17) for $`f_2`$ and (29) for $`f_1`$. Neglecting nonlinear terms in the tensors $`𝖢_i𝖯_i^{}𝟙`$, one gets

$$\zeta _1=\frac{m_1n_2}{dT_1}\sigma _{12}^{d1}v_0^3\beta _3\mu _{21}(1+\alpha _{12})\theta ^{d/2}I_\zeta (\theta ),$$
(A12)
where

$$I_\zeta (\theta )=\pi ^d𝑑𝐕_1^{}𝑑𝐕_2^{}g_{12}^{}\left[\mu _{21}g_{12}^2(1\alpha _{12})+2(𝐠_{12}^{}𝐆_{12}^{})\right]e^{\theta V_1^2V_2^2}.$$
(A13)
Here, $`𝐕_i^{}=𝐕_i/v_0`$, $`𝐠_{12}^{}=𝐠_{12}/v_0`$, $`𝐆_{12}^{}=𝐆_{12}/v_0`$, $`\theta =\mu _{12}/(\mu _{21}\gamma )`$, $`\gamma =T_1/T_2`$, and $`v_0=\sqrt{2T_2/m_2}`$. In Eqs. (A12) and (A13), use has been made of the fact that the scalar $`\zeta _1`$ cannot be coupled to the traceless tensor $`𝖢_i`$ so that the only contributions to $`\zeta _1`$ come from the (pure) Maxwellian terms in (17) and (29). The integral $`I_\zeta `$ can be evaluated by the change of variables

$$𝐱=𝐕_1^{}𝐕_2^{},𝐲=\theta 𝐕_1^{}+𝐕_2^{},$$
(A14)
with the Jacobian $`\left(1+\theta \right)^d`$. The integral $`I_\zeta (\theta )`$ can be now easily computed with the result

$$I_\zeta (\theta )=\frac{\mathrm{\Gamma }\left((d+3)/2\right)}{\mathrm{\Gamma }\left(d/2\right)}\theta ^{(d+3)/2}(1+\theta )^{1/2}\left[2\mu _{21}(1+\theta )(1+\alpha _{12})\right].$$
(A15)
Use of the result (A15) in Eq. (A12) yields

$$\zeta _1=2\frac{\pi ^{(d1)/2}}{d\mathrm{\Gamma }\left(d/2\right)}n_2\sigma _{12}^{d1}v_0\mu _{21}(1+\alpha _{12})\left(\frac{1+\theta }{\theta }\right)^{1/2}\left[2\mu _{21}(1+\theta )(1+\alpha _{12})\right].$$
(A16)
The corresponding expression for $`\zeta _2`$ can be easily obtained from (A16) and the result is

$$\zeta _2=\sqrt{2}\frac{\pi ^{(d1)/2}}{d\mathrm{\Gamma }(d/2)}n_2\sigma _2^{d1}v_0(1\alpha _{22}^2).$$
(A17)
### 2 Evaluation of $`𝖠`$ and $`𝖡`$

Since the tensor $`𝖠`$ can be easily obtained from the expression of $`𝖡`$, let us explicitly evaluate the latter tensor. It is defined as

$$𝖡=𝑑𝐕_1m_1𝐕_1𝐕_1J_{12}[𝐕_1|f_1,f_2]$$
(A18)
Using (A3), Eq. (A18) can be written as

$`𝖡`$ $`=`$ $`m_1\sigma _{12}^{d1}\beta _3\mu _{21}(1+\alpha _{12}){\displaystyle 𝑑𝐕_1𝑑𝐕_2f_1(𝐕_1)f_2(𝐕_2)g_{12}}`$ (A20)
$`\times \left[𝐠_{12}𝐆_{12}+𝐆_{12}𝐠_{12}+{\displaystyle \frac{\mu _{21}}{d+3}}(2d+33\alpha _{12})𝐠_{12}𝐠_{12}{\displaystyle \frac{\mu _{21}}{d+3}}(1+\alpha _{12})g_{12}^2𝟙\right],`$
where use has been made of (A8) and

$$𝑑\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠_{12})(\widehat{𝝈}𝐠_{12})^n\widehat{𝝈}\widehat{𝝈}=\frac{\beta _n}{n+d}g_{12}^{n2}\left(n𝐠_{12}𝐠_{12}+g_{12}^2𝟙\right).$$
(A22)
Substituting the Sonine approximations (17) and (29) for $`f_2`$ and $`f_1`$, respectively, and retaining only linear terms in the tensors $`𝖢_i`$, one gets

$`𝖡`$ $`=`$ $`m_1\sigma _{12}^{d1}\beta _3\mu _{21}(1+\alpha _{12})n_1n_2v_0^3\theta ^{d/2}\pi ^d{\displaystyle 𝑑𝐕_1^{}𝑑𝐕_2^{}g_{12}^{}e^{\theta V_1^2V_2^2}}`$ (A25)
$`\times [1+\theta 𝖢_1:(𝐕_1^{}𝐕_1^{}{\displaystyle \frac{1}{d}}V_1^2𝟙)+𝖢_2:(𝐕_2^{}𝐕_2^{}{\displaystyle \frac{1}{d}}V_2^2𝟙)]`$
$`\times \left[𝐠_{12}^{}𝐆_{12}^{}+𝐆_{12}^{}𝐠_{12}^{}+{\displaystyle \frac{\mu _{21}}{d+3}}(2d+33\alpha _{12})𝐠_{12}^{}𝐠_{12}^{}{\displaystyle \frac{\mu _{21}}{d+3}}(1+\alpha _{12})g_{12}^2𝟙\right].`$
This integral can be computed by the change of variables (A14). After a lengthy algebra, the result is

$`𝖡`$ $`=`$ $`{\displaystyle \frac{\pi ^{(d1)/2}}{d\mathrm{\Gamma }(d/2)}}\sigma _{12}^{d1}m_1n_1n_2\mu _{21}(1+\alpha _{12})v_0^3\left({\displaystyle \frac{1+\theta }{\theta }}\right)^{3/2}`$ (A28)
$`\times \{[{\displaystyle \frac{\lambda _{12}}{d+2}}+{\displaystyle \frac{d}{d+3}}\mu _{21}(1+\alpha _{12})]𝟙\mathrm{𝟚}{\displaystyle \frac{\theta }{(\mathrm{𝟙}+\theta )^\mathrm{𝟚}}}[(\mathrm{𝟙}+{\displaystyle \frac{(𝕕+\mathrm{𝟛})}{\mathrm{𝟚}(𝕕+\mathrm{𝟚})}}{\displaystyle \frac{\mathrm{𝟙}+\theta }{\theta }}\lambda _{\mathrm{𝟙𝟚}})\gamma ^\mathrm{𝟙}𝖯_\mathrm{𝟙}^{}`$
$`(1{\displaystyle \frac{(d+3)}{2(d+2)}}(1+\theta )\lambda _{12})𝖯_2^{}]\},`$
where

$$\lambda _{12}=\frac{2}{1+\theta }\frac{3}{d+3}\mu _{21}(1+\alpha _{12}).$$
(A29)
The corresponding expression for $`𝖠`$ can be easily inferred from Eq. (A28) by just making the change $`12`$ and $`\theta 1`$:

$$𝖠=\frac{\sqrt{2}\pi ^{(d1)/2}}{2d(d+2)\mathrm{\Gamma }(d/2)}\sigma _2^{d1}m_2n_2^2v_0^3(1+\alpha _{22})\left\{[d+1+\alpha _{22}(d1)]𝟙(\mathrm{𝟚}𝕕+\mathrm{𝟛}\mathrm{𝟛}\alpha _{\mathrm{𝟚𝟚}})𝖯_\mathrm{𝟚}^{}\right\}.$$
(A30)
### 3 Evaluation of $`𝛀`$

The tensor $`𝛀`$ is defined as

$$𝛀=\frac{m_1}{n_1T_1}𝑑𝐕_1𝐕_1J_{12}[𝐕_1f_{1,M},f_2].$$
(A31)
The evaluation of $`𝛀`$ can be made following similar mathematical steps as above. Thus, using (A1) and the Sonine approximation to $`f_2`$ in Eq. (A31), one gets

$`𝛀`$ $`=`$ $`{\displaystyle \frac{m_1}{T_1}}\pi ^d\sigma _{12}^{d1}\mu _{21}(1+\alpha _{12})\beta _3n_2v_0^3\theta ^{d/2}{\displaystyle 𝑑𝐕_1^{}𝑑𝐕_2^{}g^{}𝐠_{12}^{}𝐕_1^{}}`$ (A33)
$`\times [1+𝖢_2:(𝐕_2^{}𝐕_2^{}{\displaystyle \frac{1}{d}}V_2^2𝟙)]e^{\theta V_1^2V_2^2}`$
$`=`$ $`2n_2\pi ^d\mu _{21}(1+\alpha _{12})\sigma _{12}^{d1}\theta ^{d/2+1}v_0\beta _3(1+\theta )^{(1+d)}{\displaystyle 𝑑𝐱𝑑𝐲x𝐱(𝐱+𝐲)}`$ (A35)
$`\times [1+(1+\theta )^2𝖢_2:(𝐲\theta 𝐱)(𝐲\theta 𝐱)]e^{bx^2cy^2},`$
where $`b=\theta (1+\theta )^1`$ and $`c=(1+\theta )^1`$. This integral is easily performed, with the result

$$𝛀=\frac{2}{d}\frac{\pi ^{(d1)/2}}{\mathrm{\Gamma }(d/2)}\mu _{21}\sigma _{12}^{d1}n_2v_0(1+\alpha _{12})\left(\frac{1+\theta }{\theta }\right)^{1/2}\left[𝟙+\frac{\mathrm{𝟙}}{𝕕+\mathrm{𝟚}}\frac{\theta }{\mathrm{𝟙}+\theta }\left(𝖯_\mathrm{𝟚}^{}𝟙\right)\right].$$
(A36)
This expression yields Eq. (58) given in the text.