text:
# A generalization of Hall polynomials to ADE case

## 1. Hall theorem

Given a prime number $`p`$ and a positive integer $`l`$ let $`𝔽_p`$ be a field with $`p`$ elements, $`\overline{𝔽}_p`$ be an algebraic closure of $`𝔽_p`$, and $`𝔽_{p^l}`$ be the subfield of $`\overline{𝔽}_p`$ consisting of all $`x\overline{𝔽}_p`$ such that $`x^{p^l}=x`$. Thus $`𝔽_{p^l}`$ is a field with $`p^l`$ elements.

Let $`k`$ be a field, and let $`k[[t]]`$ be the ring of formal power series with coefficients in $`k`$ considered as a $`k`$-algebra. To a representation $`\pi `$ of $`k[[t]]`$ in a finite dimensional $`k`$-linear space $`D`$ one can associate a partition $`\lambda `$ as follows: $`\lambda _i=dim_k\left(\mathrm{im}\pi (t^{i1})/\mathrm{im}\pi (t^i)\right)`$. This partition is called the *type* of the representation $`\pi `$ (in \[Mac95\] the conjugate partition is called the type of $`\pi `$). Two representations are isomorphic if and only if they are of the same type.

Given a partition $`\lambda `$ let $`\lambda _1\lambda _2\mathrm{}`$ be the parts of $`\lambda `$ and let $`|\lambda |=_i\lambda _i`$. Let $`𝔒_\lambda (k)`$ be the set of all representations of $`k[[t]]`$ in $`k^{|\lambda |}`$ of type $`\lambda `$. If $`k`$ is algebraically closed then $`𝔒_\lambda (k)`$ is an irreducible smooth quasi-projective variety of dimension $`o_\lambda =_{ij}\lambda _i\lambda _j`$. Let $`𝐨_\lambda (p^l)`$ be the cardinal of $`𝔒_\lambda (𝔽_{p^l})`$. It is known (cf. \[Mac95, Chapter II\]) that $`𝐨_\lambda (p^l)`$ is a monic polynomial function of $`p^l`$ of degree $`o_\lambda `$.

Let $`\lambda `$, $`\mu ^1`$, …, $`\mu ^n`$, be partitions such that $`|\lambda |=|\mu ^1|+\mathrm{}+|\mu ^n|`$. Let $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)`$ be the set of all pairs $`(\pi ,\underset{¯}{D})`$ consisting of a representation $`\pi `$ of $`k[[t]]`$ in $`k^{|\lambda |}`$ of type $`\lambda `$ and an $`n`$-step filtration $`\underset{¯}{D}=(\{0\}=D^0D^1\mathrm{}D^n=k^{|\lambda |})`$ of $`k^{|\lambda |}`$ by subrepresentations of $`\pi `$ such that $`\pi |_{D^a/D^{a1}}`$ has type $`\mu ^a`$ for $`a=1,\mathrm{},n`$. If $`k`$ is algebraically closed then $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)`$ is a quasi-projective variety over $`k`$. The set $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (𝔽_{p^l})`$ can be identified with the set of $`𝔽_{p^l}`$-rational points in the variety $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (\overline{𝔽}_p)`$.

Steinitz \[Ste01\] and Hall \[Hal59\] observed that the sets $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)`$ have some relation to the theory of symmetric functions, and thus to the representation theory of the general linear and symmetric groups. Given a partition $`\lambda `$ such that the number of non-zero parts of $`\lambda `$ is less than $`(N+1)`$ let $`L_\lambda `$ denote the irreducible representation of $`GL(N,)`$ with the highest weight $`\lambda `$, and given an arbitrary partition $`\lambda `$ let $`\rho _\lambda `$ be the irreducible representation of $`S_{|\lambda |}`$ over $``$ associated to $`\lambda `$. Let $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ be a non-negative integer defined in either of the following two (equivalent) ways:

(1.1) $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ $`=dim_{}\mathrm{Hom}_{GL(N,)}(L_\lambda ,L_{\mu ^1}\mathrm{}L_{\mu ^n}),\text{ or}`$
(1.2) $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ $`=dim_{}\mathrm{Hom}_{S_{|\mu ^1|}\times \mathrm{}\times S_{|\mu ^n|}}(\mathrm{res}_{S_{|\lambda |}}^{S_{|\mu ^1|}\times \mathrm{}\times S_{|\mu ^n|}}\rho _\lambda ,\rho _{\mu ^1}\mathrm{}\rho _{\mu ^n}).`$
In (1.1) it is assumed that the number of non-zero parts of each of the partitions $`\lambda ,\mu ^1,\mathrm{},\mu ^n,`$ is less that $`(N+1)`$.

Let $`v^1,\mathrm{},v^n_0`$, $`u=_{a=1}^nv^a`$. Assume $`k`$ is algebraically closed. Let $`f_{v^1\mathrm{}v^n}^u=_{a<b}v^av^b`$ be the dimension of the projective variety of all $`n`$-step partial flags in $`k^u`$ with dimensions of the subfactors given by $`v^1,\mathrm{},v^n`$.

###### Theorem.

1. If $`k`$ is algebraically closed then the variety $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)`$ is empty or has pure dimension

$$dim𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)=f_{|\mu ^1|\mathrm{}|\mu ^n|}^{|\lambda |}+\frac{1}{2}\left(o_\lambda +o_{\mu ^1}+\mathrm{}+o_{\mu ^n}\right).$$
2. If $`k`$ is algebraically closed then the number of irreducible components of $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)`$ is equal to $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$.
3. There exists a computable polynomial $`𝐧_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ with integer coefficients such that the cardinal of the set $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (𝔽_{p^l})`$ is equal to $`𝐧_{\mu ^1\mathrm{}\mu ^n}^\lambda (p^l)`$ for any prime number $`p`$ and any positive integer $`l`$.
4. The polynomial $`𝐧_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ is identically equal to zero or has degree

$$\mathrm{deg}𝐧_{\mu ^1\mathrm{}\mu ^n}^\lambda =f_{|\mu ^1|\mathrm{}|\mu ^n|}^{|\lambda |}+\frac{1}{2}\left(o_\lambda +o_{\mu ^1}+\mathrm{}+o_{\mu ^n}\right).$$
5. The leading coefficient of the polynomial $`𝐧_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ is equal to $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$.

###### Remark.

$`𝐧_{\mu ^1\mathrm{}\mu ^n}^\lambda (x)=𝐨_\lambda (x)𝐡_{\mu ^1\mathrm{}\mu ^n}^\lambda (x)`$, where $`𝐡_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ is the Hall polynomial (cf. \[Mac95, Chapter II\]).

###### Proof.

Statement 1.1 is proven by Spaltenstein \[Spa82, II.5\]; 1.3, 1.4, and 1.5 are due to Steinitz \[Ste01\] and Hall \[Hal59\] (see also \[Mac95, Chapter II\]); 1.2 follows from 1.1 and 1.5. ∎

Definitions (1.1) and (1.2) of $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ are equivalent due to Schur-Weyl duality. However if one wants to generalize Theorem 1 to reductive groups other than $`GL(N)`$ there are various directions. For example certain subvarieties of parabolic flag varieties play the role of $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ for restriction multiplicities in the representation theory of Weyl groups. In this note some other varieties are considered which are relevant to tensor product multiplicities. They were independently described by Nakajima \[Nak01\], Varagnolo and Vasserot \[VV01\], and the author \[Mal01\].

## 2. A Lie algebra $`𝔤^{}`$

Let $`𝔤`$ be a simple simply laced Lie algebra over $``$ and $`𝔱`$ be a Cartan subalgebra of $`𝔤`$. Let $`𝔤^{}=𝔤𝔱`$ (a reductive Lie algebra) and let $`𝔱^{}=𝔱𝔱𝔤𝔱`$ be a Cartan subalgebra of $`𝔤^{}`$. Let $`I`$ be the set of vertices of the Dynkin graph of $`𝔤`$. Identify the weight lattice $`𝒬_𝔤`$ of $`𝔤`$ with $`[I]`$ in such a way that $`iI[I]`$ is the simple weight corresponding to the vertex $`i`$. Let $`𝒬_𝔤^{++}=_0[I]`$ be the set of positive weights. Fix a lattice $`𝒬_𝔱`$ in $`𝔱^{}`$ and identify $`𝒬_𝔱`$ with $`[I]`$ in some way. Let $`𝒬_𝔱^{++}=_0[I]`$. In what follows representations of $`𝔤^{}`$ are assumed to have the action of $`\{0\}𝔱`$ given by elements of the lattice $`𝒬_𝔱`$. One has an injective map $`\varkappa :[I][I]𝒬_𝔤𝒬_𝔱`$ given by $`\varkappa ((𝐮,𝐯))=(𝐮𝐯,𝐮𝐯+A𝐯)`$, where $`A`$ is the Cartan matrix of $`𝔤`$. Let $`𝒬_𝔤^{}`$ be the image of the map $`\varkappa `$. In what follows $`𝒬_𝔤^{}`$ is called the *weight lattice of $`𝔤^{}`$*. For $`\xi 𝒬_𝔤^{}𝒬_𝔤𝒬_𝔱`$ let $`|\xi |`$ denote the $`𝒬_𝔱`$-component of $`\xi `$. A weight $`\xi 𝒬_𝔤^{}`$ is called *integrable* if $`\varkappa ^1(\xi )_0[I]_0[I][I][I]`$. Let $`𝒬_𝔤^{}^+`$ be the set of integrable weights, and let $`𝒬_𝔤^{}^{++}=𝒬_𝔤^{}^+(𝒬_𝔤^{++}𝒬_𝔱^{++})`$. Elements of $`𝒬_𝔤^{}^{++}`$ are called *positive integrable* weights.

Given $`\lambda 𝒬_𝔤^{}^{++}`$ let $`L_\lambda `$ be an irreducible highest weight representation of $`𝔤^{}`$ with the highest weight $`\lambda `$. The representation $`L_\lambda `$ is finite dimensional and all its weights belong to $`𝒬_𝔤^{}^+`$. A finite dimensional representation of $`𝔤^{}`$ is called integrable if all its irreducible components have positive integrable highest weights. The category of finite dimensional integrable representations is closed with respect to tensor product. Let $`d_{\eta ^1\mathrm{}\eta ^n}^\xi `$ be a tensor product multiplicity:

$$d_{\eta ^1\mathrm{}\eta ^n}^\xi =dim_{}\mathrm{Hom}_𝔤^{}(L_\xi ,L_{\eta ^1}\mathrm{}L_{\eta ^n}),$$

where $`\xi ,\eta ^1,\mathrm{},\eta ^n𝒬_𝔤^{}^{++}`$ and $`|\xi |=|\eta ^1|+\mathrm{}+|\eta ^n|`$.

## 3. An algebra $`\stackrel{~}{}`$

Recall that $`I`$ is the set of vertices of the Dynkin graph of $`𝔤`$. Let $`H`$ be the set of pairs consisting of an edge of the Dynkin graph together with an orientation of the edge. Given a field $`k`$ let $`(k)`$ be the path algebra of the Dynkin graph over $`k`$, and let $``$ denote the multiplication in $`(k)`$. Fix a map $`\epsilon :H\{\pm 1\}`$ such that $`\epsilon (h)=\epsilon (\overline{h})`$, where $`\overline{h}`$ denotes the same edge as $`h`$ but with the opposite orientation. Let $`\vartheta =_{hH}\epsilon (h)h\overline{h}(k)`$, and let $`\stackrel{~}{}(k)`$ be an associative $`k`$-algebra defined as follows. As a $`k`$-linear space $`\stackrel{~}{}(k)=(k)_{iI}ku_i`$, where $`u_i`$ are some symbols, and the multiplication $``$ in $`\stackrel{~}{}(k)`$ is given by

$`f`$ $`f^{}=f\vartheta f^{},`$ $`f`$ $`u_i=f[i],`$
$`u_i`$ $`u_j=\delta _{ij}u_i,`$ $`u_i`$ $`f=[i]f,`$
where $`f,f^{}(k)`$, $`i,jI`$, $`[i](k)`$ is the path of length zero beginning and ending at a vertex $`iI`$, and $`\delta _{ij}`$ is the Kronecker symbol. The algebra $`\stackrel{~}{}(k)`$ was introduced by Lusztig. It is finitely generated as a $`k`$-algebra (cf. \[Lus00, 2.2\]).

Let $`k^I`$ be a semisimple $`k`$-algebra isomorphic to the direct sum of $`|I|`$ copies of $`k`$. In what follows $`(k)`$ and $`\stackrel{~}{}(k)`$ are considered as $`k^I`$-algebras with the embedding of the set of idempotents of $`k^I`$ into $`(k)`$ (resp. $`\stackrel{~}{}(k)`$) given by $`\{[i]\}`$ (resp. $`\{u_i\}`$).

Let $`𝐃`$ be a $`k^I`$-module (i.e. $`[I]`$-graded $`k`$-linear space). A representation of $`\stackrel{~}{}(k)`$ in $`𝐃`$ is a $`k^I`$-algebra homomorphism $`\stackrel{~}{}(k)\mathrm{End}_{k^I}𝐃`$. Let $`\dot{𝐃}=(k)_{k^I}𝐃`$ (cf. \[Lus98, Lus00\]). $`\dot{𝐃}`$ is naturally a left $`(k)`$-module. Given a representation $`\pi `$ of $`\stackrel{~}{}(k)`$ in $`𝐃`$ let $`\varpi _\pi \mathrm{Hom}_{k^I}(\dot{𝐃},𝐃)`$ be given by $`\varpi _\pi (fd)=\pi (f)d`$, and let $`𝒦_\pi `$ be the maximal $`(k)`$-submodule of $`\dot{𝐃}`$ contained in the kernel of $`\varpi _\pi `$. Let $`\xi 𝒬_𝔤^{}`$ be given by $`\xi =\varkappa ((𝐝+𝐯A𝐯,𝐯))`$, where $`𝐝=dim_{k^I}𝐃`$, $`𝐯=dim_{k^I}(\dot{𝐃}/𝒦_\pi )`$ (note that $`|\xi |=𝐝`$). Such $`\xi `$ is called the *type* of $`\pi `$. Two representations of the same type are not necessarily isomorphic. Given $`\xi 𝒬_𝔤^{}`$ let $`𝔛_\xi (k)`$ be the set of all representations of $`\stackrel{~}{}(k)`$ in $`k^{|\xi |}`$ of type $`\xi `$. It is known (cf. \[Nak98\], \[Lus98, Lus00\], \[CB00\]) that $`𝔛_\xi (k)`$ is non-empty if and only if $`\xi 𝒬_𝔤^{}^{++}`$, and that if $`k`$ is algebraically closed then $`𝔛_\xi (k)`$ is empty or is an irreducible smooth quasi-projective variety over $`k`$ of dimension $`x_\xi =_{iI}𝐯_i(2|\xi |A𝐯_i)`$, where $`(𝐮,𝐯)=\varkappa ^1(\xi )`$.

## 4. Generalized Hall theorem

Let $`\xi ,\eta ^1,\mathrm{},\eta ^n𝒬_𝔤^{}^{++}`$, $`|\xi |=|\eta ^1|+\mathrm{}+|\eta ^n|`$, and let $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)`$ be the set of all pairs $`(\pi ,\underset{¯}{𝐃})`$ consisting of a representation $`\pi `$ of $`\stackrel{~}{}`$ in $`k^{|\xi |}`$ of type $`\xi `$ and an $`n`$-step filtration $`\underset{¯}{𝐃}=(\{0\}=𝐃^0𝐃^1\mathrm{}𝐃^n=k^{|\xi |})`$ of $`k^{|\xi |}`$ by subrepresentations of $`\pi `$ such that $`\pi |_{𝐃^a/𝐃^{a1}}`$ has type $`\eta ^a`$ for $`a=1,\mathrm{},n`$. If $`k`$ is algebraically closed then $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)`$ is a quasi-projective variety over $`k`$. The variety $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)`$ is a variant of a “multiplicity variety” of \[Mal01\]. It is also implicitly contained in \[Nak01\] and \[VV01\]. The set $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (𝔽_{p^l})`$ can be identified with the set of $`𝔽_{p^l}`$-rational points in the variety $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (\overline{𝔽}_p)`$.

Let $`𝐯^1,\mathrm{},𝐯^n_0[I]`$, $`𝐮=_{a=1}^n𝐯^a`$. Assume $`k`$ is algebraically closed. Let $`g_{𝐯^1\mathrm{}𝐯^n}^𝐮=_{\begin{array}{c}a<b\\ iI\end{array}}𝐯_i^a𝐯_i^b`$ be the dimension of the projective variety of all $`n`$-step filtrations of $`k^𝐮`$ by $`k^I`$-submodules with dimensions of the subfactors given by $`𝐯^1,\mathrm{},𝐯^n`$.

###### Theorem.

1. If $`k`$ is algebraically closed then the variety $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)`$ is empty or has pure dimension

$$dim𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)=g_{|\eta ^1|\mathrm{}|\eta ^n|}^{|\xi |}+\frac{1}{2}\left(x_\xi +x_{\eta ^1}+\mathrm{}+x_{\eta ^n}\right).$$
2. If $`k`$ is algebraically closed then the number of irreducible components of $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)`$ is equal to $`d_{\eta ^1\mathrm{}\eta ^n}^\xi `$.
3. There exists a computable polynomial $`𝐩_{\eta ^1\mathrm{}\eta ^n}^\xi `$ with integer coefficients such that the cardinal of the set $`𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (𝔽_{p^l})`$ is equal to $`𝐩_{\eta ^1\mathrm{}\eta ^n}^\xi (p^l)`$ for any prime number $`p`$ and any positive integer $`l`$.
4. The polynomial $`𝐩_{\eta ^1\mathrm{}\eta ^n}^\xi `$ is identically equal to zero or has degree

$$\mathrm{deg}𝐩_{\eta ^1\mathrm{}\eta ^n}^\xi =g_{|\eta ^1|\mathrm{}|\eta ^n|}^{|\xi |}+\frac{1}{2}\left(x_\xi +x_{\eta ^1}+\mathrm{}+x_{\eta ^n}\right).$$
5. The leading coefficient of the polynomial $`𝐩_{\eta ^1\mathrm{}\eta ^n}^\xi `$ is equal to $`d_{\eta ^1\mathrm{}\eta ^n}^\xi `$.

###### Proof.

Statement 4.1 is proven in \[Nak01\] and \[Mal01\]. In these papers the base field $`k`$ is assumed to be $``$, but the proofs work for arbitrary algebraically closed field. 4.2 is proven in \[Mal01\] (it is also a corollary of \[Nak01, Lemma 4.4\] or \[VV01, Theorem 5.3\]). Statements 4.4 and 4.5 follow from 4.1, 4.2, and 4.3.

To prove 4.3 one can use the inductive argument in \[Lus00, 6.6\] replacing the variety (notation of \[Lus00\]) $`\mathrm{\Lambda }_{𝐃,𝐕}^{s,s}`$ (resp. $`\mathrm{\Lambda }_{𝐃,𝐕}^s`$) by $`𝔓_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ (resp. a tensor product variety – cf. \[Mal01\]). To prove that the number of $`𝔽_{p^l}`$-rational points of a tensor product variety depends polynomially on $`p^l`$ note (cf. \[Mal01\]) that a tensor product variety is a finite union of disjoint subsets each of which is a vector bundle over a fibration with the base equal to a $`[I]`$-graded partial flag variety and fibers isomorphic to $`\mathrm{\Lambda }_{𝐃^1,𝐕^1;𝐔^1}^s\times \mathrm{}\times \mathrm{\Lambda }_{𝐃^n,𝐕^n;𝐔^n}^s`$ for some $`𝐃^1,\mathrm{},𝐃^n`$, $`𝐕^1,\mathrm{},𝐕^n`$, $`𝐔^1,\mathrm{},𝐔^n`$ (again notation of \[Lus00\]). Thus it follows from \[Lus00, Section 6\] that the number of $`𝔽_{p^l}`$-rational points of a tensor product variety is given by a computable polynomial in $`p^l`$. ∎

###### Remark.

The Hall polynomials $`𝐡_{\mu ^1\mu ^2}^\lambda `$ (cf. Section 1) give structure constants of the Hecke algebra of $`GL(N)`$ in the basis of characteristic functions of the double cosets. The author does not know what connection (if any) exists between the polynomials $`𝐩_{\eta ^1\eta ^2}^\xi `$ and the Hecke algebra (of an extension of the algebraic group associated to $`𝔤`$ by its Cartan subgroup).

## 5. $`A_{N1}`$ case

In this section it is explained why Theorem 1 is a special case of Theorem 4. The argument is quite standard (cf. \[KP79\]). Throughout the section $`𝔤`$ is assumed to be $`𝔰𝔩_N`$, and the set of vertices $`I`$ of its Dynkin graph is identified with $`\{1,2,\mathrm{},N1\}`$ in such a way that two vertices are connected by an edge if and only if they are consequent integers.

Note that in the definition (1.1) of $`c_{\mu ^1\mathrm{}\mu ^n}^\lambda `$ the number of non-zero parts of each of the partitions $`\lambda ,\mu ^1,\mathrm{},\mu ^n,`$ is assumed to be less than $`(N+1)`$. Thus representations of $`k[[t]]`$ involved in the definition of $`𝔑_{\mu ^1\mathrm{}\mu ^n}^\lambda (k)`$ factor through $`k[[t]]/k^Nk[[t]]`$.

Recall that $`[1](k)`$ denotes the path of length zero beginning and ending at the vertex $`1`$.

###### Proposition.

Let $``$ be a two-sided ideal of $`\stackrel{~}{}(k)`$ generated by $`u_2,u_3,\mathrm{},u_{N1}`$.

1. The factor algebra $`\stackrel{~}{}(k)/`$ is generated as a $`k`$-algebra by $`(u_1+)`$ and $`([1]+)`$.
2. $`\underset{N}{\underset{}{[1][1]\mathrm{}[1]}}`$ .
3. A $`k`$-algebra homomorphism $`\phi :k[[t]]/t^Nk[[t]]\stackrel{~}{}(k)/`$ uniquely defined by $`\phi (1+t^Nk[[t]])=u_1+`$ and $`\phi (t+t^Nk[[t]])=[1]+`$ is an isomorphism.

###### Proof.

5.1 and 5.2 are straightforward. Due to 5.2 the map $`\phi `$ is well-defined, and due to 5.1 it is surjective. According to \[KP79\] there exists a faithful representation of $`k[[t]]/t^Nk[[t]]`$ which factors through $`\phi `$. Thus $`\phi `$ is injective. 5.3 follows. ∎

Note that a representation of $`\stackrel{~}{}(k)`$ factors through $`\stackrel{~}{}(k)/`$ if and only if its type $`\eta `$ satisfies the condition: $`|\eta |_i=0`$ for $`i=2,\mathrm{},N1`$. Given such $`\eta `$ let $`(𝐮,𝐯)=\varkappa ^1(\eta )`$ and let $`\nu (\eta )`$ be an $`N`$-tuple of integers given by $`\nu (\eta )_1=|\eta |_1𝐯_1`$, $`\nu (\eta )_i=𝐯_{i1}𝐯_i`$ for $`i=2,\mathrm{},N1`$, and $`\nu (\eta )_N=𝐯_{N1}`$. One has $`\eta 𝒬_𝔤^{}^{++}`$ if and only if $`\nu (\eta )`$ is an ordered partition (i.e. $`0\nu (\eta )_i\nu (\eta )_j`$ for any $`ij`$). Thus $`\nu `$ provides a bijection between types of representations of $`\stackrel{~}{}(k)/`$ and those of $`k[[t]]/t^Nk[[t]]`$. Moreover if a representation $`\pi `$ of $`\stackrel{~}{}(k)/`$ is of type $`\eta `$ then the representation $`\pi \phi `$ of $`k[[t]]/t^Nk[[t]]`$ is of type $`\nu (\eta )`$, and one has the following bijection (isomorphism of varieties if $`k`$ is algebraically closed):

$$𝔓_{\eta ^1\mathrm{}\eta ^n}^\xi (k)𝔑_{\nu (\eta ^1)\mathrm{}\nu (\eta ^n)}^{\nu (\xi )}(k).$$

Therefore Theorem 1 follows from Theorem 4.


format: text

text:
# Elasticity theory connection rules for epitaxial interfaces

## I introduction

Continuum elasticity theory was developed in the 18th and 19th centuries—prior to the general acceptance of the atomic view of matter—to describe the mechanics of elastic solids. Modern applications of elasticity theory abound throughout science and engineering, from providing a long-wavelength description of the dynamics of crystalline lattices, to the inversion of seismological data to image the three-dimensional structure of the earth’s interior.

The fundamental degree-of-freedom in a nonpolar elastic medium is the displacement field $`𝐮(𝐫)`$, the deviation of the medium at point $`𝐫`$ from its position in mechanical equilibrium. When applied to composite media consisting of layers or regions of different materials, characterized by different elastic parameters, a question naturally arises: What boundary conditions should be imposed on the displacement field at the interfaces?

An example of such a composite system is shown schematically in Fig. 1. Alternating layers of type A and type B materials, each characterized by different elastic constants and mass densities, are separated by abrupt interfaces. Within each region the displacement field satisfies an appropriate equation of motion. For an isotropic continuum with mass density $`\rho `$, the field equation is

$$\rho _t^2𝐮=v_\mathrm{l}^2\mathbf{}(\mathbf{}𝐮)v_\mathrm{t}^2\mathbf{}\times \mathbf{}\times 𝐮,$$
(1)
where $`v_\mathrm{l}\sqrt{(\lambda +2\mu )/\rho }`$ and $`v_\mathrm{t}\sqrt{\mu /\rho }`$ are the longitudinal and transverse bulk sound velocities, determined by the Lamé coefficients $`\lambda `$ and $`\mu `$. The solution of the set of second-order equations of the form (1), or its generalization to anisotropic media, requires boundary conditions on $`𝐮`$ and $`(𝐧\mathbf{})𝐮`$, where $`𝐧`$ is a unit vector normal to the interface.

The conventional boundary conditions applied in this situation (assuming fully bonded materials) are as follows: First, the displacement field is assumed to be continuous across an interface,

$$𝐮_\mathrm{A}=𝐮_\mathrm{B}.$$
(2)
The condition (2) implies that the two materials are attached and do not separate. The second condition follows from momentum conservation and requires that the force density be continuous,

$$T_\mathrm{A}^{ij}n^j=T_\mathrm{B}^{ij}n^j.$$
(3)
Here $`T^{ij}`$ is the stress tensor, defined by the continuity equation

$$_t\mathrm{\Pi }^i+_jT^{ij}=0$$
(4)
for momentum density $`𝚷\rho _t𝐮`$, and $`𝐧`$ is the unit normal. In an isotropic elastic medium, it follows from Eq. (1) that the stress tensor is given by

$`T^{ij}`$ $`=`$ $`\lambda (\mathbf{}𝐮)\delta _{ij}2\mu u_{ij}`$ (5)
$`=`$ $`c_{ijkl}u_{kl},`$ (6)
where

$$c_{ijkl}=\lambda \delta _{ij}\delta _{kl}+\mu (\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk})$$
(7)
is the elastic tensor for a linear isotropic solid, and where

$$u_{ij}(_iu_j+_ju_i)/2$$
(8)
is the strain tensor.

The purpose of this paper is to point out that these boundary conditions, (2) and (3), while quite appropriate for the geophysical application mentioned above, are generally incorrect when applied to long-wavelength vibrational dynamics in crystals with abrupt, epitaxial interfaces. The reason is because in the latter application, elasticity theory is only an approximate long-wavelength description for the underlying microscopic lattice dynamics—which necessarily depends on the detailed atomic structure of the interface—whereas (2) and (3) make no reference to that microscopic structure. For example, the correct boundary conditions must depend on the effective force constants between type A and B atoms in Fig. 1, as well as between atoms of the same type.

There are numerous applications of elasticity theory to solid state systems with heterostructures, where the use of the conventional boundary conditions would lead to quantitatively incorrect results. Examples include phonons in nanostructures such as quantum dots, quantum wells, superlattices, and nanoparticles embedded in host materials. A correct use of boundary conditions might be especially important for nanometer-scale elastic media such as phononic band-gap materials. Also, the thermal resistance of a heterojunction is determined by phonon scattering at the interface and is therefore sensitive to the connection rules or S matrix.

Finally, we would like to point out a strong analogy between this work and the problem of determining the appropriate interface boundary conditions for the envelope functions in effective mass theory. In this case, effective mass theory serves as the appropriate long-wavelength approximation to the full Schrödinger equation that contains the microscopic periodic potential of the crystalline lattice, and connection rules are required to join envelope functions through an interface between crystals with different effective mass. The microscopic theory of these connection rules was first developed by Kroemer and Zhu, and our work may be regarded as an elasticity theory analog of Refs.  and .

In the next section we give a detailed derivation of the connection rules for the case of a simple one-dimensional model of an inhomogeneous crystal, a chain of harmonic oscillators with an abrupt change in mass and spring stiffness parameters, and in Section III we compare the results of using both our new connection rules and the conventional connection rules to exact results obtained by numerical diagonalization. In Section IV we relate the connection rule problem to that of calculating the $`S`$ matrix for plane-wave scattering from the interface. The problem of determining the interface boundary conditions between three-dimensional solids is discussed in Section V, and our conclusions are summarized in Section VI.

## II connection rules in one dimension

We turn now to an analysis of the one-dimensional case, where a chain of atoms with nearest-neighbor bonds are constrained to move on a line. The vibrations in this case are purely longitudinal.

An abrupt interface is introduced at position $`x_0`$. To the left of $`x_0`$ the mass of each atom is $`m_\mathrm{A}`$, and the effective spring constant of the nearest-neighbor bonds is $`k_\mathrm{A}`$; the corresponding parameters on the right side are $`m_\mathrm{B}`$ and $`k_\mathrm{B}`$. The strength of the bond connecting the type A and B atoms, which is generally different from $`k_\mathrm{A}`$ and $`k_\mathrm{B}`$, is denoted by $`k_\mathrm{J}`$. The lattice constant on both sides is equal to $`a`$. The model we consider is illustrated in Fig. 2.

According to elasticity theory, which is valid for vibrational wavelengths large compared with $`a`$, the regions to the left and right of the interface are described by the wave equations

$$(_t^2v_\mathrm{I}^2_x^2)u_\mathrm{I}=0,v_\mathrm{I}a\sqrt{k_\mathrm{I}/m_\mathrm{I}},I=A,B.$$
(9)
The elasticity theory description of the homogeneous chain is reviewed in the appendix. To proceed, the wave equations (9) must be supplemented with boundary conditions on $`u(x_0)`$ and $`u^{}(x_0)`$.

A general linear interface boundary condition may be expressed in the form

$$\left[\begin{array}{c}u(x_0)\\ u^{}(x_0)\end{array}\right]_\mathrm{B}=M\left[\begin{array}{c}u(x_0)\\ u^{}(x_0)\end{array}\right]_\mathrm{A},$$
(10)
where $`M`$ is a $`2\times 2`$ matrix. The connection rule matrix implied by the boundary conditions (2) and (3) is

$$M=\left(\begin{array}{cc}1& 0\\ 0& k_\mathrm{A}/k_\mathrm{B}\end{array}\right).$$
(11)
A common application of (11) is to an elastic string with an abrupt change in mass density, but no change in elasticity; in this case (11) reduces to the identity matrix.

It is simple to demonstrate that (11) is the only matrix consistent with conditions (2) and (3): First, continuity requires that $`M_{11}=1`$ and $`M_{12}=0.`$ To find the other elements, we note that in one dimension the $`xx`$ component of the stress tensor (5) is $`T^{xx}=\rho v^2_xu`$. The stress immediately to the left of the interface is therefore $`T_\mathrm{A}^{xx}=k_\mathrm{A}u_\mathrm{A}^{}(x_0)`$, and that to the immediate right is $`T_\mathrm{B}^{xx}=k_\mathrm{B}u_\mathrm{B}^{}(x_0)`$. Now, Eq. (10) requires that

$$k_\mathrm{B}u_\mathrm{B}^{}(x_0)=k_\mathrm{B}\left(M_{21}u_\mathrm{A}(x_0)+M_{22}u_\mathrm{A}^{}(x_0)\right),$$
(12)
which implies

$$T_\mathrm{B}^{xx}=M_{21}k_\mathrm{B}u_\mathrm{A}(x_0)+M_{22}(k_\mathrm{B}/k_\mathrm{A})T_\mathrm{A}^{xx}.$$
(13)
Therefore, the condition (3) requires that $`M_{21}=0`$ and $`M_{22}=k_\mathrm{A}/k_\mathrm{B}`$.

We now proceed with our derivation of the correct boundary condition matrix $`M`$ for the model shown in Fig. 2. The coordinates of the atoms are written as

$$x_n(t)=x_n^0+\xi _n(t),x_n^0na.$$
(14)
The equation of motion for atom $`n`$ is

$$m_n\ddot{\xi }_n=k_\mathrm{r}(\xi _{n+1}\xi _n)k_\mathrm{l}(\xi _n\xi _{n1}),$$
(15)
where $`k_\mathrm{r}`$ is the stiffness of the spring to the right of mass $`m_n`$, and $`k_\mathrm{l}`$ is that to the left. Assuming harmonic time dependence we have, for the atoms immediately to the left ($`n=1`$) and right ($`n=0`$) of the interface,

$$\omega ^2m_\mathrm{A}\xi _1=k_\mathrm{J}(\xi _0\xi _1)k_\mathrm{A}(\xi _1\xi _2)$$
(16)
and

$$\omega ^2m_\mathrm{B}\xi _0=k_\mathrm{B}(\xi _1\xi _0)k_\mathrm{J}(\xi _0\xi _1).$$
(17)
Next we introduce the displacement field $`u(x)`$ as a smooth interpolating function between the $`\xi _n`$, such that

$$u(x_n^0)=\xi _n,$$
(18)
and use the following relations,

$`\xi _2`$ $`=`$ $`u_\mathrm{A}(x_0\frac{3}{2}a)u_\mathrm{A}(x_0)\frac{3}{2}au_\mathrm{A}^{}(x_0),`$ (19)
$`\xi _1`$ $`=`$ $`u_\mathrm{A}(x_0\frac{1}{2}a)u_\mathrm{A}(x_0)\frac{1}{2}au_\mathrm{A}^{}(x_0),`$ (20)
$`\xi _0`$ $`=`$ $`u_\mathrm{B}(x_0+\frac{1}{2}a)u_\mathrm{B}(x_0)+\frac{1}{2}au_\mathrm{B}^{}(x_0),`$ (21)
$`\xi _1`$ $`=`$ $`u_\mathrm{B}(x_0+\frac{3}{2}a)u_\mathrm{B}(x_0)+\frac{3}{2}au_\mathrm{B}^{}(x_0).`$ (22)
Because the interface boundary conditions involve the displacement field and its first derivative only, second and higher-order gradients are neglected here. Furthermore, as the frequency $`\omega `$ is formally of the order of a gradient (recall the bulk dispersion relation $`\omega =v|k|`$), for consistency we also neglect the terms proportional to $`\omega ^2`$ in Eqs. (16) and (17).

The resulting coupled equations can be put in the form

$`\left(\begin{array}{cc}k_\mathrm{J}& \frac{1}{2}ak_\mathrm{J}\\ k_\mathrm{J}& a(k_\mathrm{B}\frac{1}{2}k_\mathrm{J})\end{array}\right)`$ $`\left[\begin{array}{c}u(x_0)\\ u^{}(x_0)\end{array}\right]_\mathrm{B}`$ (23)
$`=`$ $`\left(\begin{array}{cc}k_\mathrm{J}& a(k_\mathrm{A}\frac{1}{2}k_\mathrm{J})\\ k_\mathrm{J}& \frac{1}{2}ak_\mathrm{J}\end{array}\right)`$ $`\left[\begin{array}{c}u(x_0)\\ u^{}(x_0)\end{array}\right]_\mathrm{A},`$ (24)
which, upon comparison with (10), identifies

$$\left(\begin{array}{cc}k_\mathrm{J}& \frac{1}{2}ak_\mathrm{J}\\ k_\mathrm{J}& a(k_\mathrm{B}\frac{1}{2}k_\mathrm{J})\end{array}\right)^1\left(\begin{array}{cc}k_\mathrm{J}& a(k_\mathrm{A}\frac{1}{2}k_\mathrm{J})\\ k_\mathrm{J}& \frac{1}{2}ak_\mathrm{J}\end{array}\right)$$
(25)
as the connection rule matrix. Therefore we obtain, for the model shown in Fig. 2, the connection rules

$$M=\left(\begin{array}{cc}1& a[k_\mathrm{A}k_\mathrm{B}\frac{1}{2}k_\mathrm{J}(k_\mathrm{A}+k_\mathrm{B})]/k_\mathrm{J}k_\mathrm{B}\\ 0& k_\mathrm{A}/k_\mathrm{B}\end{array}\right).$$
(26)
Several remarks are in order: First, the correct connection rules clearly depend on the microscopic structure of the interface, including the stiffness $`k_\mathrm{J}`$ of the interface bond, which is generally different than $`k_\mathrm{A}`$ and $`k_\mathrm{B}`$. The boundary conditions cannot be deduced by conservation laws that do not make reference to the microscopic structure. Second, the matrix (26) is generally off-diagonal, implying a connection between the displacement field $`u`$ on one side of the interface, with the strain $`u^{}`$, as well as the displacement, on the other. Third, the displacement field is generally not continuous through the interface, in contrast with the conventional assumption. This discontinuity, however, does not imply that the two sides are separated. It simply means that the atomic displacements $`\xi _n`$, when extrapolated from each side to the “mathematical interface” at $`x_0`$, do not meet. Fourth, we note that in the limit $`a0`$ the boundary conditions (11) and (26) agree. However, this limit is not meaningful in a real crystal. And finally, we note that (11) and (26) also become equivalent in the event that $`k_\mathrm{J}`$ has the special value $`k_\mathrm{J}^{}`$ given by

$$\frac{1}{k_\mathrm{J}^{}}=\frac{1}{2}\left(\frac{1}{k_\mathrm{A}}+\frac{1}{k_\mathrm{B}}\right).$$
(27)
## III numerical studies

When $`k_\mathrm{J}`$ differs from $`k_\mathrm{J}^{}`$, the influence of the off-diagonal element in (26) can become substantial. To demonstrate this we use elasticity theory with (11) and (26) to predict the normal modes frequencies of a one-dimensional inhomogeneous crystal of finite length $`L`$, and compare both with the exact spectrum obtained numerically. The interface is placed at $`x_0=L/2`$.

The elasticity theory spectrum is obtained by (numerically) searching for frequencies such that the three conditions $`u(0)=0`$, $`u(L)=0`$, and (10) are satisfied. The appropriate solution of the wave equation to the left of the interface, on the interval $`0xx_0`$, is

$$u_\mathrm{A}(x)=\mathrm{sin}(\omega x/v_\mathrm{A}),$$
(28)
and to the right ($`x_0xL`$) is

$$u_\mathrm{B}(x)=\alpha \mathrm{cos}(\omega x/v_\mathrm{B})+\beta \mathrm{sin}(\omega x/v_\mathrm{B}).$$
(29)
$`\alpha `$ and $`\beta `$ are uniquely determined (at each frequency) by the requirement that (10) is satisfied. This leads to

$$\left[\begin{array}{c}u(x_0)\\ u^{}(x_0)\end{array}\right]_\mathrm{B}=C\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]=M\left[\begin{array}{c}u(x_0)\\ u^{}(x_0)\end{array}\right]_\mathrm{A},$$
(30)
where

$$C\left(\begin{array}{cc}\mathrm{cos}(\omega L/2v_\mathrm{B})& \mathrm{sin}(\omega L/2v_\mathrm{B})\\ (\omega /v_\mathrm{B})\mathrm{sin}(\omega L/2v_\mathrm{B})& (\omega /v_\mathrm{B})\mathrm{cos}(\omega L/2v_\mathrm{B})\end{array}\right).$$
(31)
From (30) we obtain $`\alpha (\omega )`$ and $`\beta (\omega )`$ as

$$\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]=C^1M\left[\begin{array}{c}\mathrm{sin}(\omega L/2v_\mathrm{A})\\ (\omega /v_\mathrm{A})\mathrm{cos}(\omega L/2v_\mathrm{A})\end{array}\right],$$
(32)
and the normal mode frequencies from the remaining boundary condition $`u_\mathrm{B}(L)=0`$.

The exact spectrum is obtained by expressing the coupled equations of motion (15) for a chain of $`N`$ atoms, with the first and last atoms held fixed, as a nonsymmetric eigenvalue problem. The system size is then given by $`L=Na`$. For the results presented below, $`N=101.`$

Representative results are shown in Figs. 3 through 5. In each case the frequency $`\omega `$ of mode $`n`$ is given in units of $`\pi v_\mathrm{A}/L`$. Figs. 3 and 4 the show vibrational spectra of two inhomogeneous chains, both with $`k_\mathrm{B}=5.0k_\mathrm{A}.`$ The curves in these figures are independent of the masses $`m_\mathrm{A}`$ and $`m_\mathrm{B}`$; the only mass dependence is in the energy scale $`\pi v_\mathrm{A}/L.`$ In each case the solid line is the exact spectrum, the dotted line is the elasticity theory spectrum calculated with the conventional connection rules (11), and the dashed line is the elasticity theory spectrum calculated with the connection rules (26). In Fig. 3, $`k_\mathrm{J}=0.20k_\mathrm{A}`$, and the three spectra are similar. In Fig. 4, where $`k_\mathrm{J}=0.05k_\mathrm{A}`$, the two sides are only weakly bonded together, and the spectrum calculated with Eq. (26) agrees with the exact spectrum, whereas the spectrum calculated with Eq. (11) does not. At higher frequencies both elasticity theory spectra deviate from the exact spectrum because the wavelengths become shorter.

The final set of spectra we present, shown in Fig. 5, corresponds to a homogeneous chain, $`k_\mathrm{B}=k_\mathrm{A},`$ with a weakly bonded interface, $`k_\mathrm{J}=0.20k_\mathrm{A}`$. The spectrum calculated with Eq. (26) agrees well with the exact spectrum. The elasticity theory spectrum calculated with Eq. (11) misses the fine structure present in the exact spectrum because Eq. (11) makes no reference to the value of $`k_\mathrm{J}`$.

## IV S matrix

An alternative but physically equivalent way of expressing the interface boundary conditions is through an S matrix. Whereas the matrix $`M`$ gives the linear relation between the displacement field $`u(x_0)`$ and its derivative $`u^{}(x_0)`$ on side A to that on side B, the S matrix relates the amplitudes of waves incident on the interface, from both sides, to the corresponding outgoing waves. In this case we take $`x_0`$ to be at the origin and we write the elasticity theory solutions as

$$u_\mathrm{A}(x)=A_+e^{i\omega x/v_\mathrm{A}}+A_{}e^{i\omega x/v_\mathrm{A}}$$
(33)
and

$$u_\mathrm{B}(x)=B_+e^{i\omega x/v_\mathrm{B}}+B_{}e^{i\omega x/v_\mathrm{B}},$$
(34)
where $`A_\pm `$ and $`B_\pm `$ are complex coefficients giving the amplitudes of the plane waves shown in Fig. 6.

The S matrix relates the coefficients in Eqs. (33) and (34), and is defined by

$$\left[\begin{array}{c}A_{}\\ B_+\end{array}\right]=S\left[\begin{array}{c}A_+\\ B_{}\end{array}\right].$$
(35)
From (10) we obtain

$$\left[\begin{array}{c}B_+\\ B_{}\end{array}\right]=\left[\begin{array}{c}A_+\\ A_{}\end{array}\right]$$
(36)
and therefore

$$S=\frac{1}{_{22}}\left(\begin{array}{cc}_{21}& 1\\ \mathrm{Det}& _{12}\end{array}\right),$$
(37)
where

$$\left(\begin{array}{cc}1& 1\\ i\omega /v_\mathrm{B}& i\omega /v_\mathrm{B}\end{array}\right)^1M\left(\begin{array}{cc}1& 1\\ i\omega /v_\mathrm{A}& i\omega /v_\mathrm{A}\end{array}\right).$$
(38)
Here $`\mathrm{Det}`$ is the determinant of $``$. A useful expression for $``$ may be obtained by combining Eqs. (11) and (26) as

$$M=\left(\begin{array}{cc}1& M_{12}\\ 0& k_\mathrm{A}/k_\mathrm{B}\end{array}\right),$$
(39)
where $`M_{12}`$ is either equal to zero or to the off-diagonal element in (26). Using this representation for $`M`$ we find

$$=\frac{1}{2}\left(\begin{array}{cc}1+\frac{k_\mathrm{A}v_\mathrm{B}}{k_\mathrm{B}v_\mathrm{A}}+iM_{12}\frac{\omega }{v_\mathrm{A}}& 1\frac{k_\mathrm{A}v_\mathrm{B}}{k_\mathrm{B}v_\mathrm{A}}iM_{12}\frac{\omega }{v_\mathrm{A}}\\ 1\frac{k_\mathrm{A}v_\mathrm{B}}{k_\mathrm{B}v_\mathrm{A}}+iM_{12}\frac{\omega }{v_\mathrm{A}}& 1+\frac{k_\mathrm{A}v_\mathrm{B}}{k_\mathrm{B}v_\mathrm{A}}iM_{12}\frac{\omega }{v_\mathrm{A}}\end{array}\right)$$
(40)
and

$$\mathrm{Det}=k_\mathrm{A}v_\mathrm{B}/k_\mathrm{B}v_\mathrm{A}.$$
(41)
Note that the complex terms in the S matrix come from the off-diagonal element in (26).

The S matrix provides a simple and direct way to obtain transmission and reflection amplitudes, $`𝗍`$ and $`𝗋`$, for scattering from the interface. From (37) we observe that the transmission and reflection amplitudes for a wave of unit amplitude incident from the left $`(A_+=1\mathrm{and}B_{}=0)`$ are

$$𝗍=\frac{\mathrm{Det}}{_{22}}=\frac{2k_\mathrm{A}v_\mathrm{B}}{k_\mathrm{A}v_\mathrm{B}+k_\mathrm{B}v_\mathrm{A}iM_{12}\omega k_\mathrm{B}}$$
(42)
and

$$𝗋=\frac{_{21}}{_{22}}=\frac{k_\mathrm{A}v_\mathrm{B}k_\mathrm{B}v_\mathrm{A}iM_{12}\omega k_\mathrm{B}}{k_\mathrm{A}v_\mathrm{B}+k_\mathrm{B}v_\mathrm{A}iM_{12}\omega k_\mathrm{B}}.$$
(43)
In the limit $`k_\mathrm{A}=k_\mathrm{B}=k_\mathrm{J}`$ where the mass density is discontinuous but the elasticity is continuous, these amplitudes reduce to

$$𝗍\frac{2v_\mathrm{B}}{v_\mathrm{B}+v_\mathrm{A}}\mathrm{and}𝗋\frac{v_\mathrm{B}v_\mathrm{A}}{v_\mathrm{B}+v_\mathrm{A}},$$
(44)
the well-known results for scattering from a mass discontinuity. It can be shown that the transmission and reflection coefficients, $`𝖳`$ and $`𝖱`$, defined as the fraction of transmitted and reflected energy flux, are determined from Eqs. (42) and (43) according to

$$𝖳=\frac{v_\mathrm{A}k_\mathrm{B}}{v_\mathrm{B}k_\mathrm{A}}\left|𝗍\right|^2\mathrm{and}𝖱=\left|𝗋\right|^2.$$
(45)
In addition to relating the connection rule matrix $`M`$ to observable quantities, this scattering theory formulation serves to reemphasize the main thesis of this paper, that the connection rules must depend on the microscopic structure of the heterojunction and cannot be determined by “far field” information alone.

## V beyond one dimension

In this section we give a brief discussion of the generalization of our method to three-dimensional epitaxial heterojunctions. To allow for both longitudinal and transverse elastic waves one must work with a $`6\times 6`$ connection matrix $`\stackrel{~}{M}`$ satisfying

$$\left[\begin{array}{c}u_x(x_0)\\ u_y(x_0)\\ u_z(x_0)\\ u_x^{}(x_0)\\ u_y^{}(x_0)\\ u_z^{}(x_0)\end{array}\right]_\mathrm{B}=\stackrel{~}{M}\left[\begin{array}{c}u_x(x_0)\\ u_y(x_0)\\ u_z(x_0)\\ u_x^{}(x_0)\\ u_y^{}(x_0)\\ u_z^{}(x_0)\end{array}\right]_\mathrm{A}.$$
(46)
Here $`u_i^{}𝐧\mathbf{}u_i`$, with $`𝐧`$ a unit vector normal to the interface, and $`i=x,y,z.`$ The procedure for obtaining $`\stackrel{~}{M}`$ is identical to that described in Section II; however, in general it will be necessary to include atomic bonds beyond those connecting nearest-neighbor atoms.

To obtain quantitatively accurate connection rules one would need to determine the atomic structure of the particular interface and the required force constants. This can be accomplished using first-principles electronic structure calculation methods (for example, those based on density functional theory), although a full treatment of a three-dimensional heterojunction would be very demanding computationally.

## VI discussion

We have shown that the conventional interface boundary conditions used in elasticity theory, requiring that the displacement field and its associated stress field be continuous, are generally incorrect for epitaxial interfaces. The correct boundary conditions are nonuniversal and depend on the detailed microscopic structure of the heterojunction.

The conventional boundary conditions are incorrect because the displacement field $`𝐮(𝐫)`$ is generally discontinuous. However, this discontinuity does not imply that the two sides separate. In the elasticity theory description of crystalline lattice dynamics,

$$𝐮(𝐫_0)=𝐫_n𝐫_n^0$$
(47)
is simply a function giving the displacement of atom $`n`$ at each equilibrium lattice point $`𝐫_n^0`$. A discontinuity in $`𝐮(𝐫)`$ at a “mathematical” interface between layers of atoms implies that the atomic displacements $`𝐫_n𝐫_n^0`$ on each side of an interface do not meet when smoothly interpolated to that interface. In contrast, the condition that the stress be continuous follows from momentum conservation and is generally correct.

It is tempting to approach the interface boundary condition problem by using elasticity equations generalized to the case of a compositionally graded crystal, characterized by a position-dependent mass density and elastic parameters, and then take the limit of an abrupt composition change. But this too is incorrect, for elasticity theory is intrinsically a long-wavelength description and can be formulated only for slowly graded systems, making the required limit invalid.

For example, the generalized wave equation describing the long-wavelength vibrational dynamics in a one-dimensional crystal with lattice constant $`a`$, mass density $`\rho (x)`$, and stiffness $`k(x)`$, can be shown to be (see appendix)

$$[\rho (x)_t^2a_xk(x)_x]u(x,t)=0.$$
(48)
Integration of (48) shows that $`u(x)`$ and $`k(x)u^{}(x)`$ are continuous, consistent with the conventional boundary conditions (11). However, Eq. (48), which neglects stiffness gradients higher order than $`k^{}(x)`$, is not valid in the abrupt limit.

Having made the case that the conventional interface boundary conditions (2) and (3) do not apply to epitaxial interfaces, we must emphasize again that we have not provided generally applicable conditions to replace (2) and (3). The connection rules (26) are only valid for the simple one-dimensional interface model shown in Fig. 2.

In closing, we would like to speculate about the reason the subject of this paper has been, to the best of our knowledge, overlooked in the solid state physics literature. Historically, elasticity theory was developed as a self-contained branch of mechanics that made no reference to a possible underlying atomic structure, and much of the theory was developed before the wide acceptance of the atomic view of matter. The conventional boundary conditions (2) and (3) are certainly correct within elasticity theory proper. However, within solid state physics, elasticity theory is regarded as a long-wavelength description with a well-defined but limited regime of validity, and we believe that the connection rules in question were applied to heterostructures without considering that regime of validity.

## ACKNOWLEDGMENTS

This work was supported by the Research Corporation. It is a pleasure to thank Steve Lewis and Kelly Patton for useful discussions.

## A homogeneous chain

Here we record the long-wavelength theory of the homogeneous harmonic oscillator chain with masses $`m`$, spring constants $`k`$, and lattice constant $`a`$. In this case the equation of motion leads to

$$_t^2u(x,t)\frac{k}{m}\left[u(x+a,t)2u(x,t)+u(xa,t)\right]=0.$$
(A1)
Taylor expanding (A1) leads to the one-dimensional wave equation

$$(_t^2v^2_x^2)u(x,t)=0,$$
(A2)
with sound velocity

$$va\sqrt{k/m}.$$
(A3)
Next we derive the momentum conservation condition satisfied by the displacement field $`u`$. The momentum density carried by a longitudinal elastic wave in one-dimension is $`\mathrm{\Pi }=\rho _tu,`$ where $`\rho `$ is the mass density. In the absence of external forces, Eq. (A2) shows that $`\mathrm{\Pi }`$ satisfies the continuity equation

$$_t\mathrm{\Pi }+_xT=0,$$
(A4)
where

$$T=\rho v^2_xu$$
(A5)
is the scalar stress. As expected, (A5) is identical to the $`xx`$ component of the stress tensor (5). Similarly, the energy density $`=\frac{1}{2}\rho [(_tu)^2+v^2(_xu)^2]`$ satisfies the continuity equation

$$_t+_xj_\mathrm{e}=0,$$
(A6)
where

$$j_\mathrm{e}=\rho v^2_xu_tu$$
(A7)
is the energy flux.

The long-wavelength description of a harmonic oscillator chain with spatially varying masses and spring constants follows from the appropriate gradient expansion of

$`m(x)_t^2u(x,t)`$ $`=`$ $`k(x+\frac{a}{2})[u(x+a)u(x)]`$ (A8)
$``$ $`k(x\frac{a}{2})[u(x)u(xa)].`$ (A9)
Neglecting gradients beyond $`k^{}(x)`$ leads to the form (48) quoted in SectionVI.


text:
# Pattern formation in a metastable, gradient-driven sandpile

## I Introduction

Sandpile models were first introduced by Bak, Tang, and Wiesenfeld as explicit models of self-organized criticality (SOC) Bak et al. (1987). Since then a vast literature has analyzed sandpile properties resulting from various definitions of the toppling rules; see Refs.Mehta et al. (1996); Jensen (1998); Dhar (1999) for recent reviews.

In this paper we present the relaxation and stationary properties of a gradient-driven sandpile model with a metastable toppling rule. Interest in this type of work originates from seismology, where quasiperiodic behavior of seismic activity has been observed for certain faults and investigated with SOC-related models Dahmen et al. (1998). The model we present can describe characteristics of systems undergoing rheologic flow, e.g., in the mining environment or tectonic plates. In such systems stress can accumulate in various parts of the system for long periods of time, in contrast to normal fluid flow in which the local relaxation time is much smaller than the hydrodynamic time scale.

Furthermore, we believe that the properties we describe are of wider interest to other fields of nonequilibrium statistical mechanics where, e.g., stripe patterns similar to those that we observe appear in a variety of extended systems, including sand and biological systemsBall (1999); Shinbrot and Muzzio (2001). Our approach is quite general since it utilizes only a consistent local ordering of topplings, following the instantaneous maximum gradient, and the notion of a metastable site. Furthermore, the model we propose shows a quasiperiodic time behavior, a feature already mentioned for SOC-related models in Socolar et al. (1993) and explored recently in a similar context by ChapmanChapman (2000).

The additional feature of the proposed model is the emergence of a spatial structure for the metastable configuration driven by the transport of grains through the system. A study of this spatial configuration in and close to the stationary state is the main objective of the following sections.

The paper is organized as follows. In Sec. II we describe the toppling rule and its connection with related models. In Sec. III we present the 1D variant of the model and in Sec. IV the results for the 2D version, followed by conclusions in Sec. V.

## II The toppling rule

We consider a gradient-driven sandpile with a toppling rule which takes into account not only whether a local threshold gradient is exceeded, but also whether this situation is the result of the addition of a grain to the site under evaluation. We introduce this rule as a simple description of the dynamical weakening introduced in Ref.Dahmen et al. (1998).

The general description of the dynamics is as follows. Grains are dropped on the sites in a designated region of the lattice called the source zone; following a relaxation rule, to be described in detail subsequently, grains can change their position on the lattice until they reach an open boundary and leave the system (lattice). In short, the model describes the transport of grains from the source zone to the open boundary.

Let us now analyze in detail the toppling rule in 2D on a square lattice. At a given moment of time each site of the lattice has an associated height $`h`$ of the grain column, and we also associate with it the set of gradients $`G=\{g_l,g_r,g_u,g_d\}`$, where $`g_\alpha =hh_\alpha `$ and $`h_\alpha \{h_l,h_r,h_u,h_d\}`$ is the the height of the sand column at the four nearest neighbor sites left, right, up, and down.

We use two thresholds in our algorithm. (i)$`g_{max}`$ is the stability threshold for the maximum of $`G`$. If at least one of the gradients of $`G`$ is larger than $`g_{max}`$ the site is called metastable. It topples only if this state is the result of the receipt of a grain, either from a neighboring site which had toppled, or when a grain was dropped on the site at the start of an updating run. Metastable states do not topple when a gradient larger than $`g_{max}`$ develops as a result of the loss of grains on neighbouring sites. (ii) $`g_{min}`$ is the minimum positive gradient, which fixes the condition to stop the toppling; we call it the activity threshold. Once a site starts toppling, it sends grains, one at a time, along the instantaneous maximum gradient of $`G`$ to its nearest neighbors, provided that max$`\{G\}>g_{min}`$. If there is more than one instantaneous maximum gradient a random choice is made among them. If one site topples it will send grains to some of its nearest neighbors. We refer to the process of relaxation of one site as a ”toppling” and to the sites that have received sand grains during toppling as ”updated sites”. Their coordinates are kept in a list for the next step of the dynamics.

In one time step all the updated sites produced at the previous time step are checked for stability and, if unstable, they relax according to the above algorithm. We specify that once a site is unstable and chosen to topple the algorithm will finish the toppling sequence at the site and then will move to another site. Physically, this is equivalent with neglecting the toppling time.

Because we always topple along the instantaneous maximum gradient we can introduce a time order of the topplings. On physical grounds it seems reasonable to assume that the site that first receives a grain will topple first at the next time step, if it is unstable. We model this using a first in, first out list: the coordinates of the updated sites are stored sequentially in a list. In one time step the algorithm reads the list sequentially from the first entry: if the current site is metastable, it is toppled and the sites that are updated in this process are stored in the list for the next time step; if the current site is stable, it is simply discarded from the list. We refer to this time ordered toppling rule by the acronym TOTR. (Recall that a gradient toppling rule implies a non-Abelian sandpile where the order of toppling has definite consequences; see, e.g., Refs.Jensen (1998); Dhar (1999).)

Alternatively, we can disregard the time ordering and select at random a site from the list of updated sites for inspection of its stability. In Sec. III we make a comparative study of these two variants of the toppling rules in the 2D case. We refer to this randomly ordered toppling rule by the acronym ROTR.

In other words, we can consider one time step of the dynamics, $`\mathrm{\Delta }t`$, divided into $`N`$ bins ($`N`$ being the number of sites of the lattice), each bin containing at most the coordinates of one updated site since in a time interval $`\mathrm{\Delta }t/N`$ we can assume that typically at most one site is active. An active site generates a set of updated sites when topples. In this description the TOTR fills the bins of the next time step in the order in which the updated sites are produced at the current time step; meanwhile the ROTR fills the bins of the next time step randomly. Evidently, most of the bins are empty as the number of updated sites is much smaller than the system size and the algorithm discards the empty bins using the list of updated sites.

Thus, to resume, at a given instant of time we have three kinds of sites in our system: (i) inactive sites with all associated gradients less than $`g_{max}`$; (ii) metastable sites, which have at least one associated gradient larger than $`g_{max}`$; and (iii) active or unstable sites, which are the toppling sites at the given instant (they are perturbed metastable sites).

One observation about the times scales in the model is now in order. We can think in terms of the existence of three time scales. The smallest one is the time in which a site has toppled, $`\tau _t`$. The second time scale is the surviving time of an updated unstable site, $`\tau _s`$. The dynamics of our model is such that $`\tau _s\tau _t`$. The third time scale is the time between two droppings of grains on the lattice $`\tau _a`$. After one grain is dropped the system relaxes globally through an avalanche of topplings. Here we are ensured that $`\tau _aN_L\tau _s`$, where $`N_L`$ is the average number of time steps for the system to relax for a given size $`L`$.

We define the size of an avalanche as the total number of topplings generated in a given relaxation process after one dropping.

We stress that our toppling rule allows and introduces sites with unstable gradients after one avalanche has taken place. They emerge as the neighbors of the toppling (active) sites toward which the gradient is negative. As an unstable site topples, the negative gradient increases in absolute value, but the sites along such a direction do not receive any grains and hence are not toppled, according to our rule. When topplings in such an avalanche stop, such a site can accordingly have a maximum gradient that is larger than the threshold value $`g_{max}`$. We call such a site metastable — it can sustain gradients larger than the threshold as long as it remains unperturbed. Physically, we associate this rule with a certain local metastability of the medium through which transport takes place.

We make the observation that the metastable sites may be important in the regions where there is no dropping of grains, since they can live much longer than the average time between two droppings on the same site.

As elaborated in the next two sections, this model has properties similar to the extended version of the forest fire model described in Ref.Socolar et al. (1993), that is, quasiperiodic behavior in time and a spiked avalanche distribution.

To link further with previous studies, it is also of interest to understand the relation between our proposed model and the condition for SOC in sandpile models. References Socolar et al. (1993); Grinstein et al. (1990) present an analysis to establish the generic conditions a model has to satisfy in order to present SOC. The authors show that in two dimensions, or larger, and for conservative dynamics, an intrinsic spatial anisotropy is required to produce SOC for models which can be treated perturbatively. The model we propose has a conservative toppling rule (the noise is also conservative in the region where no dropping takes place) and evolves in two spatial dimensions under conditions of anisotropy, but it does not present the features of SOC for the avalanche distribution.

We think that the explanation of this anomaly resides in the fact that the stationary state of this model cannot be characterised as a perturbation of an interaction-free model. As we present in full detail in Sec. IV, the stationary state of this model is characterized by the appearance of deep and narrow valleys along the direction of grain transport. These features cannot be obtained as a perturbation of a diffusionlike relaxation.

We also make the observation that in the proposed toppling rules there is no external noise term, e.g., similar to the noise term in the Langevin equation. The randomness in our model comes from the selection of the toppling order in the case of the ROTR, and from the random choice of the toppling direction in the case of a degenerate maximum gradient. Another source of randomness in these models is the ‘internal noise’, in the sense discussed in Ref. Van Kampen (1992), coming from the fluctuations of the many particle dynamics which may be present even in the case of pure deterministic microscopic equations. There is no straightforward connexion between the ‘internal noise’ and the ‘external noise’ term of a Langevin equation.

In connection with this theoretical aspect we mention that recently a model with metastable states (”sticky grains” is the term used by the authors) but with a height-driven toppling rule and dissipative dynamics was studied in Ref. Mohanty and Dhar (2002) and shown to be in the directed percolation universality class.

In the following sections we present a detailed discussion of the 1D case and a statistical analysis of the mentioned patterns of metastable states for the 2D case in the stationary state and close to the stationary state.

## III The 1D case

We start our study with a one-dimensional sandpile. We choose a lattice of dimension $`L`$ with an open boundary at $`x=L`$ and with a wall at $`x=0`$. The grains are injected randomly in the region $`x[1,w]`$, with $`1wL`$, called the source zone. We choose the stability threshold $`g_{max}=2`$ and the activity threshold $`g_{min}=1`$. Let us analyze in detail the appearance of a metastable site in 1D lattice. For simplicity we start with an initial stair configuration which has a monotonic step of max$`\{G\}=2`$ descending from $`x=1`$ to $`x=L`$. This configuration is marginally stable; if we drop a grain anywhere on the lattice an avalanche occurs. First we consider a restricted source zone where grains are only dropped at $`x=1`$. A grain dropped at $`x=1`$ will then topple until $`x=L`$ for the chosen configuration. At $`x=L`$ an extra toppling takes place since $`g_{min}=1`$, and the toppling of the disturbed state continues as long as the maximum gradient is larger than $`g_{min}`$. \[Formally, the boundary condition is that $`h(L+1)=0`$, $`h(0)=\mathrm{}`$.\] When the avalanche stops, we therefore see that the gradient between the sites $`L1`$ and $`L`$ is $`3`$, hence we have a metastable site at $`L1`$. Now if we drop a new grain at $`x=1`$, the same consideration shows that the metastable configuration moves to $`L2`$. The process continues until the metastable configuration reaches the site $`x=1`$. After that a series of avalanches reconstructs the initial configuration, adding a grain at site $`L`$, then at $`L1`$, and so on.

The characteristic time period of the system is controlled by the time in which the metastable site travels from the site $`L`$ to the site $`1`$. This kind of behavior is preserved if we use a source zone with $`w>1`$ but with the ratio $`w/L`$ small. The avalanche size distribution is characterized by a peak at the end of the distribution support (see Fig. 1), which results from the avalanches produced while the metastable site is present in the transport zone between $`w+1`$ and $`L`$. The smaller size avalanches are produced after the metastable site has reached the source zone. Figure 1 shows that for fixed ration $`w/L`$ the avalanche distribution scales with $`L`$.

We close this section with the observation that the distinction between the TOTR and ROTR is irrelevant in $`1`$D since typically there is only one updated site in the algorithm.

## IV The 2D case

### IV.1 General features

In this section we compare the behavior of the system in two dimensions with the simple and well understood behavior in 1D. We choose for study a rectangular lattice of size $`L`$ with an open boundary condition at $`x=L`$ and a wall at $`x=1`$, with a periodic boundary condition along the $`y`$ axis. A simulation starts with the initial condition specified by a uniform slope of size $`1`$ per lattice step along the $`x`$ axis, the height at $`x=L`$ being $`0`$. The initial slope along the $`y`$ axis is set to $`0`$. We choose the stability threshold $`g_{max}=4`$ and the activity threshold $`g_{min}=1`$. (We have also performed simulations starting with an empty lattice; the stationary state is the same and only the transient regime is longer.) We use both toppling rules: (i) with the time ordered toppling rule (TOTR) and (ii) with a random ordered toppling rule (ROTR).

The grains are dropped randomly at $`x=1`$, $`y[1,L]`$. This source configuration allows us to explore the behavior of the transport region of the sandpile under minimal perturbation, since only one grain topples at the boundary of the transport region and the metastable sites do not have their average lifetime constrained by the average time between two droppings on the same site; hence the structure of the transport zone is influenced only by a minimal flow of grains.

The initial condition we start with can be viewed as a collection of interacting 1D sandpiles oriented along the $`x`$ axis, and consequently one might expect to see two-dimensional avalanches created by the interaction of different metastable sites.

Surprisingly we found that in the stationary state the sandpile develops a structure of valleys along the $`x`$ direction separated by terraces of sites in the metastable state. As one can see from Fig. 2 the valleys are not purely $`1`$D; instead they fluctuate slightly along the transverse direction and also show branched structures. Visual inspection of Fig. 2 shows that the transverse fluctuations and branching are more pronounced in the case of the TOTR algorithm.

After the system reaches the stationary state we notice that in the case of the TOTR the valleys do not change in time except for small fluctuations which appear toward the open boundary. In the case of the ROTR the valleys do change in time even in the stationary state. We illustrate this in Fig. 2 where the snapshots (c) and (f) follow after $`10^6`$ droppings on the snapshots (b) and (e) in the stationary regime. We observe that in the case of the TOTR the two configurations are almost identical, while for the ROTR there is a clear difference between the two configurations around the points denoted with the symbol $``$. We return to this point when we analyze the cluster size distribution and the correlation length behavior.

Another observation which we can make from visual inspection of the stationary patterns, is that the structure of the valleys changes with system size. From Fig. 3 we see that for $`L=128`$ the valleys are predominantly onedimensional. When the lattice size increases, transverse fluctuation and branches appear.

We obtain further insight about the nature of the stationary state if we analyze the time evolution of the height profile along $`z`$ (perpendicular to the $`xy`$ plane) at the boundary of the source zone. In Fig. 4 we show the time evolution of the height profile along $`y`$ for $`x=1`$ and $`x=2`$. We see that in the transient state, for both toppling rules, the source zone ($`x=1`$) has heights close to the heights of the $`x=2`$ profile. In this configuration a grain dropped in the source zone will move more probably along the $`x`$ axis, leaving the source zone in one or a few steps. In contrast, in the stationary regime we see that a completely different configuration arises. The average height of the source zone is significantly smaller than the average height of the $`x=2`$ profile. A grain dropped in the source zone will therefore first travel along the $`y`$ direction, until it meets a minimum which is connected with a valley, see Figs. 4 (b), and (e). In the stationary state the profile at the source remain unchanged for both the TOTR and ROTR as Figs. 4 (b), (c), (e), and (f) show. The data are taken from the same configuration presented in Fig. 2.

We close the presentation of general features of the model with an analysis of the avalanche size distribution in the stationary state. Figure 5 shows for both dynamics that the size distribution is not power-law-like, but rather similar to the $`1`$D case (compare Fig. 1), although with a richer structure of peaks. We remark that a scaling proportional with lattice size $`L`$ holds approximately for the plateau region and the first peak in the distribution, as in the $`1`$D case. Since the valley length scales with $`L`$, we think that this feature of the avalanche size distribution is determined by the propagation of the avalanches along the valleys.

### IV.2 Cluster characterization

Having presented a visual description of the valleys in the previous subsection, we now consider a quantitative approach. For a numerical description of the patterns shown in Figs. 2 and 3 we concentrate on the clusters formed by the bottom of the valleys. We define a site to be at the bottom of a valley if it has associated with it two or more negative gradients larger in absolute value than the threshold, i.e., at least two neighbors are metastable.

A set of bottom of the valley sites is said to form a cluster if they can be spanned by a path stepping only to nearest neighbors along the $`x`$ or $`y`$ directions.

We analyze the following quantities: the total number of clusters $`N_c(t)`$ as a function of time, the total cluster mass $`M_c(t)`$ as a function of time, the longitudinal and transverse correlation lengths $`\xi _{}`$ and $`\xi _{}`$, and, finally, the cluster size distribution at a given time, $`N(s,t)`$.

We start with the time evolution of the total number of clusters $`N_c(t)`$ and the total mass $`M_c(t)=_ssN(s,t)`$. We remind the reader that in this section a time step corresponds to a grain drop and the subsequent relaxation, and therefore we neglect the relaxation time per avalanche.

The data were collected over a single run starting with the initial condition specified at the beginning of this section. We have collected the data for one point at time intervals $`5\times 10^4`$, $`10^5`$, $`10^6`$, and $`5\times 10^6`$ for the system sizes $`L=128`$, $`256`$, $`512`$, $`1024`$, respectively. To eliminate the short time fluctuation we averaged at each point over a window of size $`10^4`$ steps from which we selected $`100`$ moments equally spaced. We mention here that we are constrained to use single runs since the characteristic time in which the system reaches the stationary state scales approximately with the third power of the lattice linear size.

Figure 6 shows that the total mass of the clusters and the number of clusters have a jump before the stationary state is reached. The jump is more pronounced in the case of the TOTR and is clearer for large system size. We observe that the time to reach the stationary state scales with the third power of system size. This fact can be explained if we assume that the average slope of the system converges to a nonzero value as the system size is increased; then the volume of the accumulated grains scales like $`L^3`$, which determines the minimal number of time steps necessary to drive the system to the stationary state.

In studies of this kind of model it is customary to test the data for scaling behaviour $`f(t,L)=L^\alpha \stackrel{~}{f}(t/L^z)`$. Our plot shows that one can define a dynamical exponent $`z=3`$ for the relaxation time, but a simple proportionality with the system size at a given power is not sustainable, although for $`L=512,1024`$ the assumption holds acceptably (the corresponding lines are parallel within a good approximation). However, the scaling hypothesis does not hold for the system size $`L=128`$, suggesting that boundary effects have a characteristic length of this order of magnitude. From Fig. 3 we see that close to the the source zone the valleys do not present branching. We note also that close to the open boundary valleys do not form, which makes it plausible that the scaling is strongly affected by finite size effects.

For a geometric charaterization of the clusters we define the following correlation lengths, following Ref. Stauffer and Aharony (1992):

$$\xi _i^2=\frac{2_sR_{i,s}^2s^2n_s}{_ss^2n_s},$$
(1)
where $`R_{i,s}^2`$ is the mean square distance along the longitudinal, $`x`$, or transverse, $`y`$, direction for the sites belonging to one cluster averaged over all clusters of size $`s`$, which can be expressed in the formula

$$2R_{i,s}^2=\frac{1}{n_s}\underset{\alpha =1}{\overset{n_s}{}}\frac{1}{s^2}\underset{m_\alpha ,n_\alpha }{}|x_{m_\alpha }^{(i)}x_{n_\alpha }^{(i)}|^2.$$
(2)
Here the index $`i=,`$ indicates the longitudinal and transverse directions, $`x_{n_\alpha }^{(i)}`$ denotes the coordinate of the sites belonging to the $`\alpha `$th cluster of size $`s`$, and $`n_s`$ is the number of clusters with size $`s`$. In other words, for each cluster from a given configuration of the lattice we compute its average square displacement from the center of mass in transverse and longitudinal directions with Eq. (2). Next we average over all the clusters from a configuration with weights specified by Eq. (1).

Figure 7 shows that the longitudinal correlation lengths jump before the stationary regime for system sizes $`L=512`$ and $`1024`$, similar to the behavior of the cluster number and the total mass( see Fig. 6).

It is intriguing that the longitudinal correlation length shows chaotic behavior in the transient state as the system size increases, especially in the case of the ROTR, even though we have averaged out the short time fluctuations. On the other hand, the time evolution of the total number of clusters does not present fluctuations of the same relative magnitude (see Fig. 6). We conjecture that the large fluctuation of the correlation lengths in the case of the ROTR is related to a high frequency of coalescence and breaking up of clusters before the stationary state is reached.

The transverse correlation presents the same kind of chaotic behavior. In any event, we see that $`\xi _{}1`$ (which is physically significant) only for the TOTR with $`L512`$. Visual inspection of Figs. 2, and 3 shows that it is this situation in which a significant number of branchings are observable.

To elucidate this point further, we study, for both types of dynamics TOTR and ROTR, the cluster size distribution at various moments of time for the lattice size $`L=1024`$. In Fig. 8 we plot the cluster distribution at two consecutive moments in the transient regime ($`t_1=8.7\times 10^8`$, $`t_2=8.75\times 10^8`$) corresponding to a minimum and a maximum of $`\xi _{}`$ observed in Fig. 7, that is, $`\xi _{}(t_1)<\xi _{}(t_2)`$, and in the stationary state ($`t_3=1.84\times 10^9`$, $`t_4=1.845\times 10^9`$) with the same condition $`\xi _{}(t_3)<\xi _{}(t_4)`$ .

The plots show that in the transient regime the states with larger $`\xi _{}`$ have a significant population of large clusters for both dynamics. In the stationary state in the case of the TOTR the cluster distribution changes very little, while for the ROTR the the large size tail presents observable fluctuations. This observation seems to point to a signature for a stationary state of the ROTR in which large scale fluctuations appear a the cluster of valleys, as shown in Fig. 2. (We have also checked this fact for other configurations.)

To summarize the two-dimensional behavior, we have studied the geometric and relaxation properties of the transport region of the sandpile. We found that as the system approaches the stationary state, clusters of deep valleys appear. Specifically, we explored the time evolution for the total number of clusters, the total mass of the clusters, the correlation lengths (longitudinal and transverse) of the clusters, and the cluster size distribution at various moments of time. The main distinction between the TOTR and ROTR is the fact that the configurations obtained from the ROTR are strongly affected by fluctuations of the longitudinal correlation length, which we have shown to be associated with fluctuation in the tail of the cluster size distribution.

The characteristic time to reach the stationary state scales with the system size as $`L^z`$, with $`z3`$, but the magnitude of the observed quantities does not follow a clear power law type of scaling. The range of data we have does not allow us to decide if this is due to strong correction from boundary effects, or if scaling is genuinely broken. The setting in of a chaotic regime at large $`L`$ complicates the analysis even further.

## V Conclusions

We have analyzed a gradient-driven sandpile model with local metastable states using two variants of the toppling rule: one which keeps the time order of the list of updated sites and one which selects a site from the list at random. We found that in 2D metastable sites generate a rough landscape with deep valleys along which grains are transported.

The valleys are organized in clusters with a pronounced development along the direction of grain flow ($`\xi _{}\xi _{}`$). For the time ordered toppling rule, we found that transverse correlations develop as the system size increases. Also, for large system size we observed chaotic behavior for the correlation lengths, associated with fluctuations in the size distribution of the clusters.

The avalanches produced in the resulting valleys have one-dimensional properties and show non-SOC behavior. We believe that this is an example of a model in the strong coupling regime which escapes the classification made in Ref. Grinstein et al. (1990) for models that can be treated perturbatively. Yet we do not have an exact mapping between our discrete model and the continuum Langevin equations used in Ref. Grinstein et al. (1990). This may be another source of the discrepancy. At the moment this is a purely numerical study of a descriptive nature. To clarify the previous uncertainties further analytical insight seems to be required. It would, e.g., be interesting to deduce a simple mean field kind of equation capable of explaining the irregular time behavior at large $`L`$ for the quantities described in the paper.

## Acknowledgments

The authors thank an anonymous referee for detailed comments and useful suggestions. One of the authors, L A, thankfully acknowledges the hospitality of the Department of Physics, Inha University, South Korea, where part of the work was done; L. A. also acknowledges partial support from the Marie Curie Fellowship of the European Community programme under Contract No. HPMF-CT-2002-01910 . We acknowledge early discussions with Morne Pistorius and preliminary calculations by him which stimulated the present investigation.


text:
# Currents and charges for mixed fields

## I Introduction

Despite being studied since long time, particle mixing and oscillations  are still among the most fascinating subjects of Particle Physics: there is a renewed interest in the physics of mixed mesons (like Kaons or B-mesons) due to the ongoing or planned experiments  and neutrino oscillations seem finally to have been observed  after a long search. On the other hand, on the theoretical side, many questions remain unanswered like the ones on the origin itself of the mixing among generations or about the values of the lepton mixing angles. Also, various issues related to the phenomenology of flavor oscillations are object of active discussions (see Refs. for recent reviews).

Recently, there has been a remarkable progress in the understanding of the Quantum Field Theory (QFT) of mixed particles, for both fermions and bosons. It has emerged  that the vacuum state for the flavor fields (describing the mixed particles) is not the same as the one for the mass eigenstates, rather is a generalized coherent state . This fact turns out to have observable consequences: the exact oscillation formulas obtained in QFT for neutrinos  and for bosons  exhibit corrections with respect to the usual textbook expressions derived in Quantum Mechanics (QM) .

The aim of this paper is to present an analysis of the currents and the charges for a system of mixed fields, which can then be used for deriving oscillation formulas. We consider both spin zero charged bosons and Dirac fermions in the simplest case of mixing among two generations.

The fact that such a problem has remained basically unsolved till now is due to the fundamental difficulties one encounters in the definition of the Hilbert space for particles with undefined mass, as the mixed particles are. Thus the results of Refs. are crucial in the following discussion. To our knowledge, the only existing treatment of currents and charges for mixed particles is the one of Ref.: the results there obtained are valid in the framework of non-relativistic Quantum Mechanics, and therefore have to be considered only approximate (i.e. they should be regarded as the non-relativistic limit of exact ones). This is so because mixing between states of different masses is not even allowed in non-relativistic QM due to superselection rules .

Apart from the intrinsic theoretical interest of our analysis, its necessity has been recently emphasized in Ref.: an unambiguous definition of currents and charges for mixed particles would lead to a more consistent way to treat particle oscillations, especially for the case (as it happens for actual experiments) in which the time of measurement is not known and what is really measured is the distance from the source. In Ref., oscillation formulas for the case of Kaon mixing are obtained after integrating the probability current density over the time of measurement and the surface of the detector ($`L`$ denotes the distance from the source and $`\alpha ,\beta `$ are flavor indices):

$`P_{\alpha \beta }(L)={\displaystyle _{t_1}^{t_2}}𝑑t{\displaystyle _A}\mathrm{𝐝𝐒}𝐣_\beta (x,t),`$ (1)
from which results in line with the usual formulas follow.

In QFT the situation is different: rather than to probability current, one has to look at expectation values of operators describing the flavor currents and charges. These expectation values have to be properly defined on the Hilbert space for mixed particles which cannot be identified with the space for the mass eigenstates, as shown in Refs.. What we show here is that currents and charges for a system of two mixed fields can be indeed properly defined in QFT: these are the currents associated with $`SU(2)`$ transformations acting on the field doublet. By using the flavor Hilbert space , we then derive exact QFT oscillation formulas (both for fermions and bosons) which agree with the ones already presented in Refs.. We do not discuss here oscillation formulas in space, although our results naturally provide a basis for such calculation, which will be treated elsewhere.

## II Currents and charges for mixed fields

Let us begin our discussion with the fermion case and consider the following Lagrangian density describing two Dirac fields with a mixed mass term (in the following we will refer explicitly to neutrinos):

$`(x)=\overline{\mathrm{\Psi }}_f(x)\left(i\overline{)}M\right)\mathrm{\Psi }_f(x),`$ (2)
where $`\mathrm{\Psi }_f^T=(\nu _e,\nu _\mu )`$ and $`M=\left(\begin{array}{cc}m_e& m_{e\mu }\\ m_{e\mu }& m_\mu \end{array}\right)`$. This Lagrangian has been already used in Ref. to model neutrino oscillations in vacuo.

The mixing transformations

$`\mathrm{\Psi }_f(x)=\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right)\mathrm{\Psi }_m(x),`$ (5)
with $`\theta `$ being the mixing angle and $`\mathrm{\Psi }_m^T=(\nu _1,\nu _2)`$, diagonalize the quadratic form of Eq.(2) to the Lagrangian for two free Dirac fields, with masses $`m_1`$ and $`m_2`$:

$`(x)=\overline{\mathrm{\Psi }}_m(x)\left(i\overline{)}M_d\right)\mathrm{\Psi }_m(x),`$ (6)
where $`M_d=diag(m_1,m_2)`$. One also has $`m_e=m_1\mathrm{cos}^2\theta +m_2\mathrm{sin}^2\theta `$, $`m_\mu =m_1\mathrm{sin}^2\theta +m_2\mathrm{cos}^2\theta `$, $`m_{e\mu }=(m_2m_1)\mathrm{sin}\theta \mathrm{cos}\theta `$.

The free fields $`\nu _i`$ (i=1,2) can be quantized in the usual way (we use $`tx_0`$):

$`\nu _i(x)={\displaystyle \underset{r}{}}{\displaystyle d^3k\left[u_{𝐤,i}^r\alpha _{𝐤,i}^r(t)+v_{𝐤,i}^r\beta _{𝐤,i}^r(t)\right]e^{i𝐤𝐱}},`$ (7)
with $`\alpha _{𝐤,i}^r(t)=e^{i\omega _{k,i}t}\alpha _{𝐤,i}^r(0)`$, $`\beta _{𝐤,i}^r(t)=e^{i\omega _{k,i}t}\beta _{𝐤,i}^r(0)`$ and $`\omega _{k,i}=\sqrt{𝐤^2+m_i^2}`$. The vacuum for the mass eigenstates is denoted by $`|0_{1,2}`$: $`\alpha _{𝐤,i}^r|0_{1,2}=\beta _{𝐤,i}^r|0_{1,2}=0`$. The anticommutation relations are the usual ones; the wave function orthonormality and completeness relations are those of Ref..

We now study the transformations acting on the doublet of free fields with different masses. $``$ is invariant under global $`U(1)`$ phase transformations of the type $`\mathrm{\Psi }_m^{}=e^{i\alpha }\mathrm{\Psi }_m`$: as a result, we have the conservation of the Noether charge $`Q=d^3xI^0(x)`$ (with $`I^\mu (x)=\overline{\mathrm{\Psi }}_m(x)\gamma ^\mu \mathrm{\Psi }_m(x)`$) which is indeed the total charge of the system (i.e. the total lepton number).

Consider now the $`SU(2)`$ transformations acting on $`\mathrm{\Psi }_m`$:

$`\mathrm{\Psi }_m^{}(x)=e^{i\alpha _j\tau _j}\mathrm{\Psi }_m(x),j=1,2,3.`$ (8)
with $`\alpha _j`$ real constants, $`\tau _j=\sigma _j/2`$ and $`\sigma _j`$ being the Pauli matrices.

For $`m_1m_2`$, the Lagrangian is not generally invariant under (8) and we obtain, by use of the equations of motion,

$`\delta (x)=i\alpha _j\overline{\mathrm{\Psi }}_m(x)[\tau _j,M_d]\mathrm{\Psi }_m(x)=\alpha _j_\mu J_{m,j}^\mu (x)`$ (9)
$`J_{m,j}^\mu (x)=\overline{\mathrm{\Psi }}_m(x)\gamma ^\mu \tau _j\mathrm{\Psi }_m(x),j=1,2,3.`$ (10)
The charges $`Q_{m,j}(t)d^3xJ_{m,j}^0(x)`$, satisfy the $`su(2)`$ algebra $`[Q_{m,j}(t),Q_{m,k}(t)]=iϵ_{jkl}Q_{m,l}(t)`$. Note that the Casimir operator is proportional to the total charge: $`C_m\left[\underset{j=1}{\overset{3}{}}Q_{m,j}^2(t)\right]^{\frac{1}{2}}`$ $`=\frac{1}{2}Q`$. From (10) we also see that $`Q_{m,3}`$ is conserved as $`M_d`$ is diagonal. We can define the combinations:

$`Q_1{\displaystyle \frac{1}{2}}Q+Q_{m,3},Q_2{\displaystyle \frac{1}{2}}QQ_{m,3}`$ (11)
$`Q_i={\displaystyle \underset{r}{}}{\displaystyle }d^3k(\alpha _{𝐤,i}^r\alpha _{𝐤,i}^r\beta _{𝐤,i}^r\beta _{𝐤,i}^r),i=1,2.`$ (12)
These are nothing but the Noether charges associated with the non-interacting fields $`\nu _1`$ and $`\nu _2`$: in the absence of mixing, they are the flavor charges, separately conserved for each generation.

Observe now that the transformation:

$`\mathrm{\Psi }_f(x)`$ $`=`$ $`e^{2i\theta Q_{m,2}(t)}\mathrm{\Psi }_m(x)e^{2i\theta Q_{m,2}(t)}`$ (13)
is just the mixing Eq.(5). Thus $`2Q_{m,2}(t)`$ is the generator of the mixing transformations: its properties have been studied in ref., where it has been shown that the action of $`G_\theta (t)e^{2i\theta Q_{m,2}(t)}`$ on the vacuum state $`|0_{1,2}`$ results in a new vector (flavor vacuum) $`|0(t)_{e,\mu }G_\theta ^1(t)|0_{1,2}`$, orthogonal to $`|0_{1,2}`$ in the infinite volume limit.

By use of $`G_\theta (t)`$, the flavor fields can be expanded as:

$`\nu _\sigma (x)`$ $`=`$ $`{\displaystyle \underset{r}{}}{\displaystyle d^3k\left[u_{𝐤,i}^r\alpha _{𝐤,\sigma }^r(t)+v_{𝐤,i}^r\beta _{𝐤,\sigma }^r(t)\right]e^{i𝐤𝐱}}.`$ (14)
with $`(\sigma ,i)=(e,1),(\mu ,2)`$. The flavor annihilation operators are defined as $`\alpha _{𝐤,\sigma }^r(t)G_\theta ^1(t)\alpha _{𝐤,i}^r(t)G_\theta (t)`$ and $`\beta _{𝐤,\sigma }^r(t)G_\theta ^1(t)\beta _{𝐤,i}^r(t)G_\theta (t)`$.

Let us now perform the $`SU(2)`$ transformations on the flavor doublet $`\mathrm{\Psi }_f`$:

$`\mathrm{\Psi }_f^{}(x)=e^{i\alpha _j\tau _j}\mathrm{\Psi }_f(x),j=\mathrm{\hspace{0.17em}1},2,3.`$ (15)
We obtain

$`\delta (x)=i\alpha _j\overline{\mathrm{\Psi }}_f(x)[\tau _j,M]\mathrm{\Psi }_f(x)=\alpha _j_\mu J_{f,j}^\mu (x),`$ (16)
$`J_{f,j}^\mu (x)=\overline{\mathrm{\Psi }}_f(x)\gamma ^\mu \tau _j\mathrm{\Psi }_f(x),j=\mathrm{\hspace{0.17em}1},2,3.`$ (17)
The related charges $`Q_{f,j}(t)`$ $``$ $`d^3xJ_{f,j}^0(x)`$, still close the $`su(2)`$ algebra and $`C_f=C_m=\frac{1}{2}Q`$. Due to the off–diagonal (mixing) terms in the mass matrix $`M`$, $`Q_{f,3}(t)`$ is time–dependent. This implies an exchange of charge between $`\nu _e`$ and $`\nu _\mu `$, resulting in the flavor oscillations.

In accordance with Eq.(11), we define the flavor charges for mixed fields as

$`Q_e(t)`$ $``$ $`{\displaystyle \frac{1}{2}}Q+Q_{f,3}(t)`$ (18)
$`Q_\mu (t)`$ $``$ $`{\displaystyle \frac{1}{2}}QQ_{f,3}(t)`$ (19)
with $`Q_e(t)+Q_\mu (t)=Q`$.

These charges have a simple expression in terms of the flavor ladder operators ($`\sigma =e,\mu `$):

$`Q_\sigma (t)`$ $`=`$ $`{\displaystyle \underset{r}{}}{\displaystyle d^3k\left(\alpha _{𝐤,\sigma }^r(t)\alpha _{𝐤,\sigma }^r(t)\beta _{𝐤,\sigma }^r(t)\beta _{𝐤,\sigma }^r(t)\right)},`$ (20)
which happens because they are connected to the Noether charges $`Q_i`$ of Eq.(11) via the mixing generator: $`Q_\sigma (t)=G_\theta ^1(t)Q_iG_\theta (t)`$. Notice also that the operator

$`\mathrm{\Delta }Q_e`$ $``$ $`Q_e(t)Q_1=Q_{f,3}(t)Q_{m,3}`$ (21)
describes how much the mixing violates the lepton number conservation for a given generation.

The oscillation formulas are obtained by taking expectation values of the above charges on the (flavor) neutrino state. Consider for example an initial electron neutrino state defined as $`|\nu _e\alpha _{𝐤,e}^r(0)|0_{e,\mu }`$ (for a discussion on the correct definition of flavor states see Refs.). Working in the Heisenberg picture, we obtain

$`{}_{e,\mu }{}^{}0|Q_\sigma (t)|0_{e,\mu }^{}=\mathrm{\hspace{0.17em}0},`$ (22)
$`𝒬_{𝐤,\sigma }(t)\nu _e|Q_\sigma (t)|\nu _e`$ (23)
$`=\left|\{\alpha _{𝐤,\sigma }^r(t),\alpha _{𝐤,e}^r(0)\}\right|^2+\left|\{\beta _{𝐤,\sigma }^r(t),\alpha _{𝐤,e}^r(0)\}\right|^2.`$ (24)
where $`|0_{e,\mu }|0(0)_{e,\mu }`$. Charge conservation is obviously ensured at any time: $`𝒬_{𝐤,e}(t)+𝒬_{𝐤,\mu }(t)=\mathrm{\hspace{0.33em}1}`$. We remark that the expectation value of $`Q_\sigma `$ cannot be taken on vectors of the Fock space built on $`|0_{1,2}`$, as shown in Refs..

The oscillation (in time) formulas for the flavor charges, on an electron neutrino state, then follow:

$`𝒬_{𝐤,e}(t)=\mathrm{\hspace{0.17em}1}`$ $``$ $`\mathrm{sin}^2(2\theta )|U_𝐤|^2\mathrm{sin}^2\left({\displaystyle \frac{\omega _{k,2}\omega _{k,1}}{2}}t\right)`$ (25)
$``$ $`\mathrm{sin}^2(2\theta )|V_𝐤|^2\mathrm{sin}^2\left({\displaystyle \frac{\omega _{k,2}+\omega _{k,1}}{2}}t\right),`$ (26)
with $`|U_𝐤|^2+|V_𝐤|^2=1`$ and

$`|V_𝐤|={\displaystyle \frac{|𝐤|}{2\sqrt{\omega _{k,1}\omega _{k,2}}}}\left(\sqrt{{\displaystyle \frac{\omega _{k,1}+m_1}{\omega _{k,2}+m_2}}}\sqrt{{\displaystyle \frac{\omega _{k,2}+m_2}{\omega _{k,1}+m_1}}}\right)`$ (27)
This result is exact. The differences with respect to the usual formula for neutrino oscillations are in the energy dependence of the amplitudes and in the additional oscillating term. For $`|𝐤|\sqrt{m_1m_2}`$, we have $`|U_𝐤|^21`$ and $`|V_𝐤|^20`$ and the traditional (Pontecorvo) oscillation formula is recovered.

Finally, it is useful to define the flavor currents in the following way:

$`J_e^\alpha (x)`$ $``$ $`{\displaystyle \frac{1}{2}}I^\alpha (x)+J_{f,3}(x),`$ (28)
$`J_\mu ^\alpha (x)`$ $``$ $`{\displaystyle \frac{1}{2}}I^\alpha (x)J_{f,3}(x).`$ (29)
The sum of these two currents gives $`I^\alpha `$, i.e. the Noether current density for the total system. These quantities are the natural candidates for a QFT formulation of flavor oscillations in space.

For the bosonic case we can proceed in a similar way. For simplicity, we deal with spin zero stable particles. The Lagrangian for two mixed charged scalar fields is (where necessary, we use a “hat” symbol to distinguish similar quantities from the fermionic case):

$`(x)`$ $`=`$ $`_\mu \mathrm{\Phi }_f^{}(x)^\mu \mathrm{\Phi }_f(x)\mathrm{\Phi }_f^{}(x)\widehat{M}\mathrm{\Phi }_f(x)`$ (30)
$`=`$ $`_\mu \mathrm{\Phi }_m^{}(x)^\mu \mathrm{\Phi }_m(x)\mathrm{\Phi }_m^{}(x)\widehat{M}_d\mathrm{\Phi }_m(x)`$ (31)
with $`\mathrm{\Phi }_f^T=(\varphi _A,\varphi _B)`$ being the flavor fields and $`\widehat{M}=\left(\begin{array}{cc}m_A^2& m_{AB}^2\\ m_{AB}^2& m_B^2\end{array}\right)`$. Those are connected to the free fields $`\mathrm{\Phi }_m^T=(\varphi _1,\varphi _2)`$ with $`\widehat{M}_d=diag(m_1^2,m_2^2)`$ by a rotation similar to Eq.(5). For an extensive treatment of boson mixing in QFT, see Ref..

Similarly to the fermionic case, we consider the transformation $`\mathrm{\Phi }_m^{}(x)=e^{i\alpha _j\tau _j}\mathrm{\Phi }_m(x)`$ with $`j=1,2,3`$, which gives

$`\delta (x)=i\alpha _j\mathrm{\Phi }_m^{}(x)[\tau _j,\widehat{M}_d]\mathrm{\Phi }_m(x)=\alpha _j_\mu \widehat{J}_j^\mu (x)`$ (32)
$`\widehat{J}_{m,j}^\mu (x)=i\mathrm{\Phi }_m^{}(x)\tau _j\stackrel{}{^\mu }\mathrm{\Phi }_m(x),j=\mathrm{\hspace{0.17em}1},2,3,`$ (33)
where $`\stackrel{}{^\mu }\stackrel{}{^\mu }\stackrel{}{^\mu }`$.

The charges $`\widehat{Q}_{m,i}(t)`$ $`d^3x\widehat{J}_{m,j}^0(x)`$, satisfy the $`su(2)`$ algebra (at each time $`t`$). $`2\widehat{Q}_{m,2}(t)`$ is the generator of the boson mixing transformations whose properties are studied in Refs.. The “flavor” vacuum is defined as $`|0(t)_{_{A,B}}\widehat{G}_\theta ^1(t)|\widehat{0}_{1,2}`$ with $`\widehat{G}_\theta (t)e^{i2\theta \widehat{Q}_{m,2}(t)}`$ and the flavor annihilation operators are $`a_{𝐤,A}(t)\widehat{G}_\theta ^1(t)a_{𝐤,1}\widehat{G}_\theta (t)`$, etc.. As in the fermionic case, these operators involve Bogoliubov coefficients, which now satisfy the hyperbolic relation $`|\widehat{U}_𝐤|^2|\widehat{V}_𝐤|^2=1`$ with

$`|\widehat{V}_𝐤|{\displaystyle \frac{1}{2}}\left(\sqrt{{\displaystyle \frac{\omega _{k,1}}{\omega _{k,2}}}}\sqrt{{\displaystyle \frac{\omega _{k,2}}{\omega _{k,1}}}}\right).`$ (34)
The flavor charges for bosons are:

$`\widehat{Q}_\sigma (t)`$ $`=`$ $`{\displaystyle d^3k\left(a_{𝐤,\sigma }^{}(t)a_{𝐤,\sigma }(t)b_{𝐤,\sigma }^{}(t)b_{𝐤,\sigma }(t)\right)},`$ (35)
where $`\sigma =A,B`$. We define a flavor $`\mathrm{`}\mathrm{`}A\mathrm{"}`$ state as $`|a_Aa_{𝐤,A}^{}(0)|0_{_{A,B}}`$. On this state, the expectation values of the flavor charges give

$`\widehat{𝒬}_\sigma (t)`$ $``$ $`a_A|\widehat{Q}_\sigma (t)|a_A`$ (36)
$`=`$ $`\left|[a_{𝐤,\sigma }(t),a_{𝐤,A}^{}(0)]\right|^2\left|[b_{𝐤,\sigma }^{}(t),a_{𝐤,A}^{}(0)]\right|^2.`$ (37)
We also have $`{}_{_{A,B}}{}^{}0|\widehat{Q}_\sigma (t)|0_{_{A,B}}^{}=0`$ and $`\widehat{𝒬}_{𝐤,A}(t)+\widehat{𝒬}_{𝐤,B}(t)=1`$. The explicit calculation gives

$`\widehat{𝒬}_{𝐤,A}(t)=\mathrm{\hspace{0.17em}1}`$ $``$ $`\mathrm{sin}^2(2\theta )|\widehat{U}_𝐤|^2\mathrm{sin}^2\left({\displaystyle \frac{\omega _{k,2}\omega _{k,1}}{2}}t\right)`$ (38)
$`+`$ $`\mathrm{sin}^2(2\theta )|\widehat{V}_𝐤|^2\mathrm{sin}^2\left({\displaystyle \frac{\omega _{k,2}+\omega _{k,1}}{2}}t\right).`$ (39)
For $`|𝐤|^2\frac{m_1^2+m_2^2}{2}`$ we have $`|\widehat{V}_𝐤|0`$ and the usual result is recovered. Notice the negative sign in front of $`|\widehat{V}_𝐤|^2`$ in contrast with the one for fermions: the oscillating flavor charge for bosons can assume also negative values (of course, charge conservation holds for the overall system of two mixed fields). This happens since in QFT, flavor states are essentially multiparticle ones and one has to abandon the single–particle probabilistic interpretation of QM in favor of a statistical picture.

It is important to stress again that the use of the flavor Hilbert space, i.e. of the Hilbert space constructed on the flavor vacuum, is absolutely crucial in obtaining the above oscillation formulas Eqs.(26), (39). In fact, although the currents and the charges for the flavor fields are defined purely from symmetry considerations on the Lagrangian, the choice of the representation is essential when their expectation values have to be defined. It turns out indeed, that arbitrary mass parameters can be introduced in the expansion of the flavor fields, and that the only quantities free from such arbitrary quantities are the expectation values of the flavor charges on states defined as above in the flavor Hilbert space.

## III Conclusions

In this paper, we have shown that it is possible to define consistently currents and charges for a system of mixed fields, i.e. for fields with no definite masses. This was done in the case of two flavors, by studying the $`SU(2)`$ transformations acting on a doublet of fields with different masses: in this way, we naturally recovered the generator of the mixing transformations already introduced in Refs. for fermions and bosons, respectively. Our analysis also leads in a straightforward way to the exact QFT oscillation formulas found recently  which exhibit relevant corrections with respect to the usual textbook QM ones.

Our results provide a basis for an exact QFT calculation of oscillation formulas with dependence on space rather than time, a situation which is much closer to the experimental setup. This is one of the aspect which we intend to study in the future.

In the present paper we considered two–flavor mixing for simplicity. The above analysis can be extended to the mixing among an higher number of generations. The case of three–flavor fermion mixing is particularly relevant for the phenomenology (quarks, neutrinos) and such an extension is currently under investigation.

This work has been partially supported by INFN and MURST.