Bi-continuum modelling of layered structures and crystalline interfaces
bi-cont1.tex
\documentstyle[12pt]{article}

\topmargin=0in
\headheight=0pt
\headsep=0pt
\textwidth 6.5in      %skutkuje, implikuje te� \hsize
\textheight 9.0in      %skutkuje, implikuje te� \vsize
\oddsidemargin 0pt    %skutkuje
\def\gt{\ifmmode{{\mathchar"313E}}\else{$\rangle$}\fi}
\def\lt{\ifmmode{{\mathchar"313C}}\else{$\langle$}\fi}


\begin{document}
\title {Bi-continuum modelling of layered structures and
crystalline interfaces}
\author {Ma{\l}gorzata Sztyren\thanks{Department of Mathematics and Information Science,
Warsaw University of Technology, Pl. Politechniki 1, PL-00-661 Warsaw.
E--mail: emes@mech.pw.edu.pl}
}
\date{}
\maketitle
\hfill To the memory of Ekkehart Kr\"{o}ner\\

\input{epsf.tex}
\input{trans.tex}
\special{ps:}

\newcommand {\eqn}[1]{\begin{equation}#1\end{equation}}
\newcommand {\eqna}[1]{\begin{eqnarray}#1\end{eqnarray}}
\newcommand {\la}{\longrightarrow}

{\footnotesize {\bf Summary: }
The bi-continuum model composed of two interpenetrating and
dynamically coupled material continua is analysed as a simplified but
relatively accurate way to describe some physical phenomena
in crystalline solids. The essential novelty of our approach
consists in treating a crystalline medium as a bi-continuum,
even if the crystalline lattice is structurally
single-component. Particular attention is paid to the oscillatory
behaviour of solutions on the atomic level.
Starting from a discrete atomic chain, the basic formulation of
the bi-continuum model is derived. The essential features of the
model, including accuracy of the results
as functions of physical parameters, are discussed. 
}

\tableofcontents

\newpage
\section{Introduction}
\indent
    Discrete structural models such as atomic lattices or
stratified structures are capable of delivering, as a rule, highly
relialable and accurate results concerning a wide range of physical
phenomena in solid matter. Most frequently, however, to produce
such results, an atomic model has to be treated with the aid of
computer simulation techniques, what imposes severe restrictions
on the size of computed samples and leads to the lack of
generality in the results. Moreover, atomic models are sensitive
to details of the interatomic potentials; the
final reliability of results depends critically on the
reliability of the potentials used. In consequence,
simplified models which may produce analytical -- although
less accurate -- results are still highly desirable. The loss of
accuracy should not, however, be catastrophic.

    In the present paper we focus our attention on atomic
configurations in complex layered structures or in the
vicinity of planar or quasi-planar crystalline interfaces
\cite{MS Udine}.
In such situations, the configurations frequently exhibit
a sign-alternating behaviour 
\cite{Bristowe+Crocker,Rogula+Czerwinska 2001}
on atomic scale,
not reproducible by a simple classic continuum. We have found,
however, that such a behaviour
can be described by an appropriately chosen two-component continuum
model, called the bi-continuum.

   The concept of multi-continuum, understood as a material
medium on which several
displacement fields are defined, is intuitively evident when one
considers multicomponent media (like e.g. mixture of substances
with different physical properties
or the medium composed of electrons considered as one component and
heavy ions as the other one).

   The idea given in \cite{copenhagen}
was to introduce the multicontinuum
description for structurally complex media, even if such a medium
has only one component. The different displacement fields are
defined on different substructures. We shall show how the
idea of multicontinuum works in the case when multi=2.

   From the point of view of independent variables, the problem is
essentially one-dimensional. The corresponding variable will be
denoted by \(x\) with, according to the needs, appropriate indices
or diacritics. For the sake of simplicity, we perform our theoretical
constructions on linear atomic chains, although in the real
three-dimensional applications the "atoms" are intended to represent
whole crystalline planes or appropriate structural layers.

   In Sections 2 and 3 we state the problem for a homogeneous
atomic chain and examine the sign-alternating solutions near the
"interface". In Section 4 we analyse a composite linear chain,
consisting of two components. This provides the basis for construction
of the bi-continuum model in Section 5. Subsequently we solve
the interface problem in atomic and bi-continuum models, respectively,
compare the solutions, and discuss the resulting degree of
approximation.


\section{Homogeneous atomic chain}\label{Kajtek}
\indent
    Let us consider an atomic chain with interactions between the first and
second nearest neighbours. Let $u_{i}$  denote the displacement of the atom
number $i$ from its equilibrium position. By  $\kappa_{f}, n_{f}$
and $\kappa_{s}, n_{s}$ we denote the elastic constant and initial stretch
connected with shorter and longer bonds, respectively.\\
\vskip 0pt
\noindent \hfil\hbox{\epsffile{fig1.ps}}\hfil\\
\vskip -0.5\baselineskip
%
The energy of a short bond is given by the formula
\eqn{W_{f}=\frac{1}{2}\kappa_{f}(\nabla_{1}u)^{2}+n_{f}\nabla_{1}u}
while that of a long one by
\eqn{W_{s}=\frac{1}{2}\kappa_{s}(\nabla_{2}u)^{2}+n_{s}\nabla_{2}u,}
where $\nabla_{1}$ and $\nabla_{2}$  are symbols of the central
difference operators acting on displacements of the first and second
neighbours, respectively.
The elastic energy of the chain equals the sum over all the bonds.
\par
      The state of equilibrium of a homogeneous chain is described by
the set of equations
\eqn{\kappa_{f}(u_{i-1}-2u_{i}+u_{i+1})+\kappa_{s}(u_{i-2}-2u_{i}+u_{i+2})=0,
           \label{ll3}}
where the index $i=1,2,3,\ldots$ .
The general solution to equations (\ref{ll3}) is given by the function
\eqn{ u_{i}=Az^{i}+Bz^{-i}+Ei+C}
with arbitrary constants $A, B, E, C$,  and
\eqn{ z = -\beta+\sqrt{\beta^{2}-1},\label{l3}}
where
\eqn{ \beta=\frac{\alpha}{2}+1,\ \  \alpha=\frac{\kappa_{f}}{\kappa_{s}}.}
\section{Atomic chain with "interface"\label{sec4}}

   Now, instead of a homogeneous chain we shall consider two
semi--infinite homogeneous parts (of possibly different materials)
connected by an `interface'
modelled by intermediate bonds of the third kind. The quantities related
to the intermediate bonds will be labelled by the superscript $^{o}$, while
the superscripts $^{-}$ and $^{+}$ will refer to the left and right
semi--infinite chains, respectively.\\
\par
For simplicity we assume that the parameter $\alpha$
has the same value for all three materials,
\eqn{ \frac{\kappa^{-}_{f}}{\kappa^{-}_{s}}=
      \frac{\kappa^{o}_{f}}{\kappa^{o}_{s}}=
      \frac{\kappa^{+}_{f}}{\kappa^{+}_{s}}=\alpha .}
Moreover, we assume that\\

a) the homogeneous chains are in equilibrium,\\

b) there are no forces at infinity.\\

The above conditions imply the following relation between
the initial stretch forces:
\eqn{ n_{f}+2n_{s}=0,}
which means that the sum of forces transmitted by any section vanishes.
\par
      The solution of such an interface problem is given by
\eqn{u_{j}=\left\{ \begin{array}{ll} \displaystyle
          \frac{\varepsilon_{-}}{1-z}z^{-j-1} & {\rm for } \ j=-1,-2,\dots\\
          \ & \ \\
          \displaystyle
          A-\frac{\varepsilon_{+}}{1-z}z^{j-1} & {\rm for } \ j=1,2,\dots\ ,
            \end{array} \right .\label{l9}}
where $z$  is defined by the formula (\ref{l3}). Displacements are calculated
relative to the asymptotic value at  $-\infty$. The asymptotic value
at $+\infty$  is given by
\eqn{A=\varepsilon_{o}+\frac{\varepsilon_{-}+\varepsilon_{+}}
                            {1-z} .}
The symbols  $\varepsilon_{-},\varepsilon_{o}$ and $\varepsilon_{+}$  denote
deformations and are expressed by material parameters as follows:
\eqn{\varepsilon_{-}=-T_{1}(\varepsilon_{o}+\frac{1}{\kappa_{s}^{o}}
                           (n_{s}^{o}-n_{s}^{-}))\ \ ,}
\eqn{\varepsilon_{o}=\frac{1}{\kappa_{s}^{o}}(2\beta-T_{1}-T_{2})^{-1}
     (T_{1}(n_{s}^{o}-n_{s}^{-})-T_{2}(n_{s}^{+}-n_{s}^{o}))\ \ ,
     \label{eps0}}
\eqn{\varepsilon_{+}=T_{2}(\varepsilon_{o}-\frac{1}{\kappa_{s}^{o}}
                           (n_{s}^{+}-n_{s}^{o}))\ \ ,}
where we have introduced the auxiliary notations
\eqn{T_{1}=(1-\frac{\kappa_{s}^{-}}{\kappa_{s}^{o}}(1+\frac{1}{z}))^{-1}\ \ ,}
\eqn{T_{2}=(1-\frac{\kappa_{s}^{+}}{\kappa_{s}^{o}}(1+\frac{1}{z}))^{-1}\ \ .}
The resulting solution is represented by sharp-bend vertices of the line in
Fig. 3.

   Let us notice that the number  $z$  is always negative.
It implies a strongly oscillating character of the solution given
by formula (\ref{l9}): the nearest neighbours are displaced in opposite
directions.
Thus every second atom is displaced in the same direction.
Hence, as a result
of the presence of the second neighbour interactions, we have obtained
a polarization of the chain, expressed as a relative displacement
of both substructures.

\section{Composite atomic chain}

To construct the bi-continuum model we start from introduction of two fields
of displacements: $v$ and $w$, each of them defined on one of the sub-chains.
\vskip 7pt
\noindent \hfil\hbox{\epsffile{fig2.ps}}\hfil\\
\vskip -0.5\baselineskip
%
We shall describe by the indices \(x\) and \(\bar{x}\) the position of "black"
and "white" atoms supporting the fields \(v\) and \(w\), respectively. The
description may be applied to two-component structures as well as to
one-component but dimerized ones. Altogether, in the translationally
invariant case we
have four classes of bonds. Two bonds
between the first neighbours are charcterised by parameters \(\kappa_f,\ n_f\)
and \(\bar{\kappa}_f,\ \bar{n}_f\), while two other bonds, between the second
neighbours,
are characterised by \(\kappa_s,\ n_s\) and \(\bar{\kappa}_s,\ \bar{n}_s\).
If there are no forces at infinity,
then the sum of forces transmitted by any section vanishes,
%
\eqna{n_s+n_f+\bar{n}_s=0,\nonumber\\
      n_s+\bar{n}_f+\bar{n}_s=0.}
This fact implies, in turn, that initial stretches of bonds between
the nearest neighbours equal each other:
%
\eqn{n_f=\bar{n}_f.\label{nf}}
%
The symbol \(a\) denotes now the distance between the closest atoms of
the same colour (which corresponds to \(2a\) in the previous sections),
\(b\) and \(c\) denote the distances between atoms of different colours.
These quantities are connected by the relation
%
\eqn{b+c=a.}
The energy of the chain can be composed of 4 parts, each one connected with
one kind of bonds,
%
\eqn{W=W_s+W_{\bar{s}}+W_f+W_{\bar{f}},}
where
\eqna{W_s=\sum(\frac{1}{2}\kappa_s(v_x-v_{x-a})^2+
   n_s(v_x-v_{x-a})),\nonumber \\
   W_{\bar{s}}=\sum(\frac{1}{2}\bar{\kappa}_s(w_{\bar{x}+a}-w_{\bar{x}})^2+
   \bar{n}_s(w_{\bar{x}+a}-w_{\bar{x}}))),\nonumber \\
   W_f=\sum(\frac{1}{2}\kappa_f(v_x-w_{\bar{x}})^2+
   n_f(v_x-w_{\bar{x}})),\nonumber \\
   W_{\bar{f}}=\sum\frac{1}{2}\bar{\kappa}_f(w_{\bar{x}}-v_{x-a})^2+
   \bar{n}_f(w_{\bar{x}}-v_{x-a})).
   \label{wwww}}
%
The last two parts of the energy contain terms combining values of different
fields at different points. We shall express these terms with the help of
fields \(\bar{v}\) and \(\bar{w}\) which interpolate the fields
\(v\) and \(w\) to the whole chain. By definition,
%
\eqn{\bar{v}_{\bar{x}}=\frac{\kappa_fv_x+\bar{\kappa}_fv_{x-a}}
                {\kappa_f+\bar{\kappa}_f},\ \ \ \bar{v}_x=v_x}
and
\eqn{\bar{w}_x=\frac{\bar{\kappa}_fw_{\bar{x}+a}+\kappa_fw_{\bar{x}}}
            {\kappa_f+\bar{\kappa}_f},\ \ \ \
           \bar{w}_{\bar{x}}=w_{\bar{x}}.}
%
Using these interpolating fields, we can express the energy of the first
neighbour bonds by
%
\eqn{W_f=\sum(\frac{1}{2}\kappa_f[(v_x-\bar{w}_x)+
        \frac{\bar{\kappa}_f}{\kappa_f+\bar{\kappa}_f}
        (w_{\bar{x}+a}-w_{\bar{x}})]^2+n_f[(v_x-\bar{v}_x)+
        \frac{\bar{\kappa}_f}{\kappa_f+\bar{\kappa}_f}
        (w_{\bar{x}+a}-w_{\bar{x}})]),}
\eqn{W_{\bar{f}}=\sum(\frac{1}{2}\kappa_{\bar{f}}
      [(w_{\bar{x}}-\bar{v}_{\bar{x}})+
        \frac{\kappa_f}{\kappa_f+\bar{\kappa}_f}
        (v_x-v_{x-a})]^2+\bar{n}_f[(w_{\bar{x}}-\bar{v}_{\bar{x}})+
        \frac{\kappa_f}{\kappa_f+\bar{\kappa}_f}
        (v_x-v_{x-a})]).\label{ww}}
%
All the differences are now expressed either by two fields at the same
point, or by one field at different points.

\section{The bi-continuum}

We arrive at the corresponding bi-continuum expression for energy
by appropriate truncations of the Taylor series. As the simplest
possibility let us consider the correspondence
%
\eqna{(1)\ \ \ \nabla{f} \la af' \ ,\nonumber\\
      (2)\ \ \sum \la\ \ \frac{1}{a}\int\ ,\nonumber\\
      (3)\ \ \ \bar{f} \la f.\ \ \ \ \ \ \ \
      \label{cor0}}
The symbol $f$ represents here an arbitrary function, the dash denotes
the derivative with respect to \(x\), and $a$ -- the lattice parameter.
The set of rules
(\ref{cor0}) will be referred to as the correspondence of order 0.
We shall consider also
another set of rules, called correspondence of order 1, which takes
into account the microstructure of the elementary cell of our composite
chain. Namely, instead of rule (3) in (\ref{cor0}), we take
%
\eqna{(3a)\ \ \ \bar{v}\la\ \ v+\frac{\bar{\kappa}}
      {\kappa_f+\bar{\kappa}_f}v',\nonumber\\
      (3b)\ \ \bar{w}\la\ w-\frac{\bar{\kappa}}
      {\kappa_f+\bar{\kappa}_f}w',}
%
where 
\eqn{\bar{\kappa}=\kappa_fc-\bar{\kappa}_fb.\label{barka}}
%
If \(\bar{\kappa}_f=\kappa_f\) and
\(b=c\), then both the correspondence rules are equivalent.
Let us note that the transition to correspondence rules of order 1 does not
augment the order of equations. 
\par
    Using the correspondence rule of order 0 one obtains the following
formulae for energies of the corresponding bonds:
%
\eqn{W_s=\int( \frac{a}{2}\kappa_sv'^2+n_sv')dx,}\nonumber


\eqn{W_{\bar{s}}=\int(\frac{a}{2}\bar{\kappa}_sw'^2+\bar{n}_fw')dx,}\nonumber


\eqn{W_f=\frac{1}{2a}\int[\kappa_f
  (v-w+\frac{\bar{\kappa}_f}{\kappa_f+\bar{\kappa}_f}aw')^2+
  an_f(v-w+\frac{\bar{\kappa}_f}{\kappa_f+\bar{\kappa}_f}aw')]dx,}\nonumber

\eqn{W_{\bar{f}}=\frac{1}{2a}\int[\bar{\kappa}_f
  (w-v+\frac{\kappa_f}{\kappa_f+\bar{\kappa}_f}av')^2+
  an_f(w-v+\frac{\kappa_f}{\kappa_f+\bar{\kappa}_f}av')]dx.}
%
In the last expression the equality (\ref{nf}) has been used. The energy
connected with the n.n. bonds may be rewritten in the form
\eqna{W_f+W_{\bar{f}}=\frac{1}{2a}\int[(\kappa_f+\bar{\kappa}_f)
  (v-w)^2+a^2\frac{\kappa_f\bar{\kappa}_f}{(\kappa_f+
  \bar{\kappa}_f)^2}\nonumber
  (\bar{\kappa}_fw'^2+\kappa_fv'^2)\\
  +\frac{an_f}{\kappa_f+
  \bar{\kappa}_f}(\bar{\kappa}_fw'+\kappa_fv')]dx
  +boundary\ term.}
%
The correspondence rule of order 1 introduce definite corrections to
the above formulae. These corrections result in
\par
a) an additional coupling term of the form \((v-w)(v'+w')\),
\par
b) renormalisation of coefficients in other terms.

   In consequence, the energy density of the homogeneous bi-continuum
takes the form
\eqn{w_d=\frac{1}{2}ev'^2+\frac{1}{2}gw'^2+\bar{\kappa}(v-w)(v'+w')+
         \frac{1}{2}h(v-w)^2,}
%
with the phenomenological coefficients. One can connect them with parameters
of the discrete model by the following equalities:
%
\eqna{e=a\kappa_s+a\bar{\kappa}_f(\frac{\kappa_f-\bar{\kappa}}
       {\kappa_f+\bar{\kappa}_f})^2,\nonumber\\
      g=a\bar{\kappa}_s+a\kappa_f(\frac{\bar{\kappa}_f-\bar{\kappa}}
       {\kappa_f+\bar{\kappa}_f})^2,\nonumber\\
      h=\frac{{\kappa}_f+\bar{\kappa}_f}{a},\label{pheno}}
%
and \(\bar{\kappa}\) is given by the formula (\ref{barka}).

\par
      The equations of static equilibrium for a homogeneous medium
have the form
%
\eqna{ev''-2\bar{\kappa}w'-h(v-w)=0,\nonumber\\
      gw''+2\bar{\kappa}v'+h(v-w)=0,}
%
and, without forces at infinity, have the general solution
%
\eqna{v=\left\{ \begin{array}{ll}
        Ae^{-\zeta x}+C, & {\rm for} \ x{\gt}0,\\
        \ & \ \\
        Be^{\zeta{x}} & {\rm for } \ x{\lt}0, \nonumber\\
        \end{array} \right\} \nonumber\\
      w=\left\{ \begin{array}{ll}
         Ak_1e^{-\zeta x}+C, & {\rm for} \ x{\gt}0,\\
        \ & \ \\
        Bk_2e^{\zeta{x}} & {\rm for } \ x{\lt}0.
        \end{array} \right\}
         \label{solu}}
%
with the exponent coefficient expressed by material constants according to
the formula
%
\eqn{\zeta^2=\frac{h(e+g)-4\bar{\kappa}^2}{eg}.}
%
The coefficients connecting amplitudes of fields \(v\) and \(w\) equal
%
\eqn{k_1=\frac{h-\zeta^2e}{h-2\zeta\bar{\kappa}},\ \ \
     k_2=\frac{h-\zeta^2e}{h+2\zeta\bar{\kappa}}  .}
%
In the case of exact inversion symmetry, when \(\kappa_s=\bar{\kappa}_s, b=c, 
\kappa_f=\bar{\kappa}_f\) and then \(\bar{\kappa}=0\),
%
\eqn{\zeta^2=\frac{2h}{e},\ \ k_1=k_2=-1.\label{dz}}
%
In the next section we shall consider such a symmetric case.

\section{Comparison between bi-continuum and discrete solutions}
   The solution of the bi-continuum interface problem may be constructed
from functions
(\ref{solu}), with appropriate values of constants \(A, B, C\). 
The simplest way to match this solution to its discrete counterpart is to
request that values of corresponding branches of functions \(u, v\) and
\(w\) coincide at the atoms nearest to the interface. Two such
conditions
plus coincidence of asymptotic values will be sufficient to calculate
\(A, B\) and \(C\). In this way, making use of (\ref{eps0}), we have
%
\eqn{C=\frac{T_1(n^0_s-n^-_s)}{\kappa^0_s(\beta-T_1)},}
%
and, for the other constants
%
\eqn{A=B=-\frac{\varepsilon_1}{1-z}e^{\zeta a/4}.}

\newpage
\noindent
The bi-continuum fields \(v\) and \(w\) are plotted in Fig. 3.
\vskip 7pt
\noindent \hfil\hbox{\epsffile{fig3.ps}}\hfil\\
\vskip -0.5\baselineskip
%
\noindent
The upper branches correspond to the field \(w\).
The values of
discrete function \(u\) correspond to the sharp-bends
of the broken line.
\vskip 7pt
\noindent \hfil\hbox{\epsffile{fig4.ps}}\hfil\\
\vskip -0.5\baselineskip
%
   To examine the accuracy of the bi-continuum model,
the error in value of fields \(v\) and \(w\) at the distance \(3a/4\)
from the interface, equivalent to the second nearest neighbours,
is plotted vs. \(1/\alpha\) in Fig. 4.
The error corresponding to \(\alpha=0.05\) equals about 4 percent of
the lattice constant.
Note that, under the assumed correspondence rules, the error
does not change the sign.

   Another way to compare the atomic and bi-continuum models is to
analyse the relation between the local dimerisation ranges
or, equivalently, between the exponents in both models.
Let us note that, since the symbol \(a\) denotes the distance between
the second n.n., one should calculate the discrete exponent \(\lambda\)
from the equation
%
\eqn{z^2=e^{-\lambda a},}
%
which, according to (\ref{l3}), implies the dependence of \(\lambda\)
on \(\alpha\). On the other hand, for the bi-continuum model, by
making use of equations (\ref{dz}) and (\ref{pheno}), one obtains
%
\eqn{\zeta^2=\frac{16}{a^2}\frac{\alpha}{4+\alpha}.\label{ze}}
%
The above equation implies \(0{\lt}a\zeta{\lt}4\). 
\vskip 7pt
\noindent \hfil\hbox{\epsffile{fig5.ps}}\hfil\\
\vskip -0.5\baselineskip
%

\noindent
The dependence of \(1/\lambda\) on \(1/\zeta\), is plotted in Fig. 5
as the lower curve to be compared with the line \(\lambda=\zeta\). A
substantial difference is present only at distances smaller than
the lattice constant \(a\).
\newpage
\noindent
To enable a more precise estimation of the error, the difference
\(1/\lambda - 1/\zeta\)
is plotted vs. \(1/\alpha\) in Fig. 6.
\vskip 7pt
\noindent \hfil\hbox{\epsffile{fig6.ps}}\hfil\\
\vskip -0.5\baselineskip
%
\section{Final remarks}
   We have presented a method of construction of the bi-continuum model
based on a heuristic system - a linear atomic chain. One may, however,
forget about the derivation, and consider the obtained equations as
phenomenological ones. Generally, it is not necessary to relate the
phenomenological coefficients to the parameters of an atomic model.
The fairly good agreement between the bi-continuum and atomistic solutions
can be treated as an argument in behalf of our model, at least in the
sense that we have not made gross errors.

   The internal structure of the "atoms" in the chain may be complex --
they need not be interpreted as material points. As an important example
one can mention atomic planes of 3D crystalline structures, in particular
the oxygen-copper planes, or even the octahedral layers in high
temperature superconducing materials \cite{Rogula+Czerwinska 2001}.
   If our "atoms" are planes or layers \cite{Kroner 67},
one should replace, in the
expression for
energy, the scalar quantities by apropriate vectors and matrices.
And the construction will run analogically, under the condition that one
takes into account the interactions of the first and second n.n.
between planes or layers.
   Similarly one can consider layered structures in the form
of arrays of Josephson junctions. Such systems occur in a natural
way in superconducting crystals, \cite{Walker+Luettmer 96}
or else they can be shaped artificially
in a technological way.

   How important are the planar interfaces? There are many cases
where they are important by themselves. It is expedient to notice
that (what we call) the "planar" interfaces need not
be strictly planar. It is sufficient that the curvatures of the
corresponding surfaces should be not too high as compared with the
characteristic exponents (both of dimension \(cm^{-1}\),
inverse length).

   The interfaces can also affect the bulk properties of materials.
It concerns the materials which contain a lot of interface area per
unit volume. Also this property can be expressed by a parameter of
dimension of inverse length which is particularly high for
nanostructured materials.

\section{Acknowledgements}
This work was supported by the State Committee for Scientific Research
(Poland) under grant No. 7 TO7A 010 16.\\
\noindent
The paper has been presented at the International Symposium
on Structured Media, Pozna\a'{n}, 16-21 September 2001


\thebibliography{12}
\bibliography{}
\end{document}


epsf.tex
%   EPSF.TEX macro file:
%   Written by Tomas Rokicki of Radical Eye Software, 29 Mar 1989.
%   Revised by Don Knuth, 3 Jan 1990.
%
%   TeX macros to include an Encapsulated PostScript graphic.
%   Works by finding the bounding box comment,
%   calculating the correct scale values, and inserting a vbox
%   of the appropriate size at the current position in the TeX document.
%
%   To use, simply say
%   \input epsf           % somewhere early on in your TeX file
%   \epsfbox{filename.ps} % where you want to insert a vbox for a figure
%
%   The effect will be to typeset the figure as a TeX box, at the
%   point of your \epsfbox command. By default, the graphic will have its
%   `natural' width (namely the width of its bounding box, as described
%   in filename.ps). The TeX box will have depth zero.

%   You can enlarge or reduce the figure by saying
%     \epsfxsize=<dimen> \epsfille{filename.ps}
%   instead. Then the width of the TeX box will be \epsfxsize, and its
%   height will be scaled proportionately.
%   (The \epsfbox macro resets \epsfxsize to zero after each use.)

%   If you want TeX to report the size of the figure (as a message
%   on your terminal when it processes each figure), say `\epsfverbosetrue'.

%
\newread\epsffilein    % file to \read
\newif\ifepsffileok    % continue looking for the bounding box?
\newif\ifepsfbbfound   % success?
\newif\ifepsfverbose   % report what you're making?
\newdimen\epsfxsize    % horizontal size after scaling
\newdimen\epsfysize    % vertical size after scaling
\newdimen\epsftsize    % horizontal size before scaling
\newdimen\epsfrsize    % vertical size before scaling
\newdimen\epsftmp      % register for arithmetic manipulation
\newdimen\pspoints     % conversion factor

\pspoints=1truebp      % Adobe points are `big'
\epsfxsize=0pt         % Default value, means `use natural size'

%
\def\epsfbox#1{%
%
%   The first thing we need to do is to open the
%   PostScript file, if possible.
%
\openin\epsffilein=#1
\ifeof\epsffilein\errmessage{I couldn't open #1, will ignore it}\else
%
%   Okay, we got it. Now we'll scan lines until we find one that doesn't
%   start with %. We're looking for the bounding box comment.
%
   {\epsffileoktrue \chardef\other=12
    \def\do##1{\catcode`##1=\other}\dospecials \catcode`\ =10
    \loop
       \read\epsffilein to \epsffileline
       \ifeof\epsffilein\epsffileokfalse\else
%
%   We check to see if the first character is a % sign;
%   if not, we stop reading (unless the line was entirely blank);
%   if so, we look further and stop only if the line begins with
%   `%%BoundingBox: '.
%
          \expandafter\epsfaux\epsffileline. \\%
       \fi
   \ifepsffileok\repeat
   \ifepsfbbfound\else
      \ifepsfverbose\message{No bounding box comment in #1; using defaults}\fi
      \global\def\epsfllx{72}%
      \global\def\epsflly{72}%
      \global\def\epsfurx{540}%
      \global\def\epsfury{720}\fi
   }\closein\epsffilein
%
%   Now we have to calculate the scale and offset values to use.
%   First we compute the natural sizes.
%
   \epsfrsize=\epsfury\pspoints
   \advance\epsfrsize by-\epsflly\pspoints
   \epsftsize=\epsfurx\pspoints
   \advance\epsftsize by-\epsfllx\pspoints
%
%   If `epsfxsize' is 0, we default to the natural size of the picture.
%   Otherwise we scale the graph to be \epsfxsize wide.
%
   \ifnum\epsfxsize=0 \epsfxsize=\epsftsize \epsfysize=\epsfrsize
%
%   We have a sticky problem here:  TeX doesn't do floating point arithmetic!
%   Our goal is to compute y = rx/t. The following loop does this reasonably
%   fast, with an error of at most about 16 sp (about 1/4000 pt).
%   
   \else\epsftmp=\epsfrsize \divide\epsftmp\epsftsize
     \epsfysize=\epsfxsize \multiply\epsfysize\epsftmp   
     \multiply\epsftmp\epsftsize \advance\epsfrsize-\epsftmp
     \epsftmp=\epsfxsize
     \loop \advance\epsfrsize\epsfrsize \divide\epsftmp 2
     \ifnum\epsftmp>0
        \ifnum\epsfrsize<\epsftsize\else
           \advance\epsfrsize-\epsftsize \advance\epsfysize\epsftmp \fi
     \repeat     
   \fi
%
%  Finally, we make the vbox and stick in a \special that dvips can parse.
%
   \ifepsfverbose\message{#1: width=\the\epsfxsize, height=\the\epsfysize}\fi
   \epsftmp=10\epsfxsize \divide\epsftmp\pspoints
   \vbox to\epsfysize{\vfil\hbox to\epsfxsize{%
      \special{psfile=#1 llx=\epsfllx\space lly=\epsflly\space
          urx=\epsfurx\space ury=\epsfury\space rwi=\number\epsftmp}%
      \hfil}}%
\fi\epsfxsize=0pt}%

%
%   We still need to define the tricky \epsfaux macro. This requires
%   a couple of magic constants for comparison purposes.
%
{\catcode`\%=12 \global\let\epsfpercent=%\global\def\epsfbblit{%BoundingBox:}}%
%
%   So we're ready to check for `%BoundingBox:' and to grab the
%   values if they are found.
%
\long\def\epsfaux#1#2 #3\\{\ifx#1\epsfpercent
   \def\testit{#2}\ifx\testit\epsfbblit
      \epsfgrab #3 . . . \\%
      \epsffileokfalse
      \global\epsfbbfoundtrue
   \fi\else\ifx#1\par\else\epsffileokfalse\fi\fi}%
%
%   Here we grab the values and stuff them in the appropriate definitions.
%
\def\epsfgrab#1 #2 #3 #4 #5\\{\global\def\epsfllx{#1}\global\def\epsflly{#2}%
   \global\def\epsfurx{#3}\global\def\epsfury{#4}}%
%
%   Finally, another definition for compatibility with older macros.
%
\let\epsffile=\epsfbox

trans.tex
%---------------------------------------------------------------
% PostScript transformations for DVIPS        (c) BOP 1993, 1994
% BOP s.c.,  81-805 Sopot, al. Niepodleg�o�ci 783, tel. 51-23-09
%---------------------------------------------------------------
% USAGE:
%
%   1) 
%      x_scale and y_scale are given in percents:
%      \zscale changes uniformly x_size and y_size (x_scale=y_scale=scale);
%      \xscale changes only x_size (y_scale=100);
%      \yscale changes only y_size (x_scale=100);
%
%         \zscale{scale} followed by an \hbox or a \vbox
%         \xyscale{x_scale}{y_scale} followed by an \hbox or a \vbox
%         \yxscale{y_scale}{x_scale} followed by an \hbox or a \vbox
%         \xscale{x_scale} followed by an \hbox or a \vbox
%         \yscale{y_scale} followed by an \hbox or a \vbox
%
%   2) 
%      \xyscaleto resizes a box uniformly such that the resulting wd=x_dimen;
%      \yxscaleto resizes a box uniformly such that the resulting ht=y_dimen;
%      \xscaleto resizes a box horizontally such that the resulting wd=x_dimen;
%      \yscaleto resizes a box verically such that the resulting ht=y_dimen;
%
%         \xyscaleto{x_size} followed by an \hbox or a \vbox
%         \yxscaleto{y_size} followed by an \hbox or a \vbox
%         \xscaleto{x_size} followed by an \hbox or a \vbox
%         \yscaleto{y_size} followed by an \hbox or a \vbox
%
%   3)
%      \revolve rotates anticlockwise the box (either \hbox or \vbox
%      following the command) by 90 degree; for the resulting box
%      width = height + depth of the original box, height = width
%      of the original box, and the reference point is the left top corner
%      of the original box; this means that revolving a box four times
%      yields the original box if and only if its depth is zero.
%      \revolvedir- is equivalent to \revolve, \revolvedir+ does almost
%      the same, but rotates the box clockwise;
%     
%         \revolve followed by an \hbox or a \vbox
%         \revolvedir+ followed by an \hbox or a \vbox
%         \revolvedir- followed by an \hbox or a \vbox
%     
%   4) 
%      \rotate rotates a box bay an arbitrary angle, clockwise for angle>0,
%      the resulting box width=height=depth=0pt;
%
%         \rotate{angle} followed by an \hbox or a \vbox
%
%   5)
%      \xflip and \yflip flip the box horizontally and vertically, resp.,
%      i.e., with respect to a vertical and horizontal axis of the box,
%      without changing dimensions of the \box;
%     
%         \xflip followed by an \hbox or a \vbox
%         \yflip followed by an \hbox or a \vbox
%     
%   6) 
%      \slant slopes a box by the angle alpha such that tan(alpha)=slant,
%      without changing dimensions of the \box;
%
%         \slant{slant} followed by an \hbox or a \vbox
%     
%---------------------------------------------------------------
% HISTORY:
% 18 VIII 1993 ver. 0.1
%    * first release
% 30--31 VIII 1993 ver. 0.2
%    * third parameter eliminated from \scale (via \afterassignment
%      and \aftergroup hackery)
%    * added \zscale, \xyscale, \yxscale, \xscale, and \yscale
%      with scaling given in percents
%    * added \scaleto, \xyscaleto, \yxscaleto, \xscaleto, and \yscaleto
% 3 IX 1993 ver. 0.21
%    * \the_scale renamed to \lastscale and made global, thus it became 
%      available to a user
% 8 IX 1993 ver. 0.22
%    * all transformations return \hbox, because of currentpoint
%      positioning
%    * the old version of \scale is become undefined
%    * \slant and \rotate fit the new convention of parameter's
%      hackery -- they are assumed to be followed by an \hbox
%      or a \vbox
% 6 XII 1993 ver. 0.23
%    * \revolve added
%    * \rotate with \vbox patched
% 10 II 1994 ver. 0.24
%    * \revolvedir+, \revolvedir-, \xflip, and \yflip  added
%---------------------------------------------------------------
\edef\undtranscode{\the\catcode`\_} \catcode`\_11
%---------------------------------------------------------------
\newbox\box_tmp % temporary box register
\newdimen\dim_tmp % temporary dimen register (for arithmetic manipulation)
%---------------------------------------------------------------
\def\jump_setbox{\aftergroup\after_setbox}% a general trick
%---------------------------------------------------------------
% ``floating point arithmetic'' (excerpted from T. Rokicki):
% r  y
%
% ^
% |
% |
% |
% |
% |
% 0--------------> t   x
%
\def\resize
    % dimen registers:
    #1% y   make y such that y/r=x/t
    #2% r
    #3% x
    #4% t
%   We have a sticky problem here:  TeX doesn't do floating point arithmetic!
%   Our goal is to compute y = rx/t. The following loop does this reasonably
%   fast, with an error of at most about 16 sp (about 1/4000 pt).
    {%
    % save parameters to the internal variables:
    \dim_r#2\relax \dim_x#3\relax \dim_t#4\relax
    \dim_tmp=\dim_r \divide\dim_tmp\dim_t
    \dim_y=\dim_x \multiply\dim_y\dim_tmp
    \multiply\dim_tmp\dim_t \advance\dim_r-\dim_tmp
    \dim_tmp=\dim_x
    \loop \advance\dim_r\dim_r \divide\dim_tmp 2
    \ifnum\dim_tmp>0
      \ifnum\dim_r<\dim_t\else
        \advance\dim_r-\dim_t \advance\dim_y\dim_tmp \fi
    \repeat
    % assign result:
    #1\dim_y\relax
    }
\newdimen\dim_x    % horizontal size after scaling
\newdimen\dim_y    % vertical size after scaling
\newdimen\dim_t    % horizontal size before scaling
\newdimen\dim_r    % vertical size before scaling
%\newdimen\dim_tmp % register for arithmetic manipulation (already declared)
%---------------------------------------------------------------
\def\perc_scale#1#2{% #1 -- xscale, #2 -- yscale, in percents,
                    % to be followed by an \hbox or a \vbox
  \def\after_setbox{%
    \hbox\bgroup
    \special{ps:  gsave
             currentpoint #2 100 div div exch #1 100 div div exch
             currentpoint neg exch neg exch
             translate #1 100 div #2 100 div scale translate}%
    \dim_tmp\wd\box_tmp \divide\dim_tmp100 \wd\box_tmp#1\dim_tmp
    \dim_tmp\ht\box_tmp \divide\dim_tmp100 \ht\box_tmp#2\dim_tmp
    \dim_tmp\dp\box_tmp \divide\dim_tmp100 \dp\box_tmp#2\dim_tmp
    \box\box_tmp \special{ps: grestore}\egroup}%
    \afterassignment\jump_setbox\setbox\box_tmp =
}%
%---------------------------------------------------------------
{\catcode`\p12 \catcode`\t12 \gdef\PT_{pt}}
\def\hull_num{\expandafter\hull_num_}
\expandafter\def\expandafter\hull_num_\expandafter#\expandafter1\PT_{#1}
%---------------------------------------------------------------
\def\find_scale#1#2{% #1 -- size after rescaling, #2 -- \wd or \ht
% Finds a scale (\lastscale macro) such that the box following the macro
% call would have the respective dimen (i.e., #2) equal to #1 after rescaling.
% NOTE: it is assumed that prior to calling \find_scale a macro
%       \extra_complete is defined.
  \def\after_setbox{%
    \resize\dim_tmp{100pt}{#1}{#2\box_tmp}%
    \xdef\lastscale{\hull_num\the\dim_tmp}\extra_complete}%
  \afterassignment\jump_setbox\setbox\box_tmp =
}
%---------------------------------------------------------------
\def\scaleto#1#2#3#4{% #1 -- size of dimen #2 (\wd or \ht) after scaling
                     % #3 -- actual x-size, #4 -- actual y-size
  \def\extra_complete{\perc_scale{#3}{#4}\hbox{\box\box_tmp}}%
  \find_scale{#1}#2}
%---------------------------------------------------------------
\let\xyscale\perc_scale
\def\zscale#1{\xyscale{#1}{#1}}
\def\yxscale#1#2{\xyscale{#2}{#1}}
\def\xscale#1{\xyscale{#1}{100}}
\def\yscale#1{\xyscale{100}{#1}}
%---------------------------------------------------------------
\def\xyscaleto#1{\scaleto{#1}\wd\lastscale\lastscale}
\def\yxscaleto#1{\scaleto{#1}\ht\lastscale\lastscale}
\def\xscaleto#1{\scaleto{#1}\wd\lastscale{100}}
\def\yscaleto#1{\scaleto{#1}\ht{100}\lastscale}
%---------------------------------------------------------------
\def\slant#1{% #1 (slant) = tan(alpha), where alpha is the "italic" angle.
             % To be followed by an \hbox or a \vbox
  \hbox\bgroup
  \def\after_setbox{%
    \special{ps: gsave 0 currentpoint neg exch pop 0 currentpoint exch pop
             translate [1 0 #1 1 0 0] concat translate}%
    \box\box_tmp \special{ps: grestore}\egroup}%
  \afterassignment\jump_setbox\setbox\box_tmp =
}%
%---------------------------------------------------------------
\def\rotate#1{% #1 -- angle,
              % to be followed by an \hbox or a \vbox
  \hbox\bgroup
  \def\after_setbox{%
    \setbox\box_tmp\hbox{\box\box_tmp}% otherwise does not work with \vbox
    \wd\box_tmp 0pt \ht\box_tmp 0pt \dp\box_tmp 0pt
    \special{ps: gsave currentpoint currentpoint translate
             #1 rotate neg exch neg exch translate}%
    \box\box_tmp
    \special{ps: grestore}%
    \egroup}%
  \afterassignment\jump_setbox\setbox\box_tmp =
}%
%---------------------------------------------------------------
\newdimen\box_tmp_dim_a
\newdimen\box_tmp_dim_b
\newdimen\box_tmp_dim_c
\def\plus_{+}
\def\minus_{-}
\def\revolvedir#1{% to be followed by an \hbox or a \vbox
  \hbox\bgroup
% check parameter:
   \def\param_{#1}%
   \ifx\param_\plus_ \else \ifx\param_\minus_
   \else
     \errhelp{I would rather suggest to stop immediately.}%
     \errmessage{Argument to \noexpand\revolvedir should be either + or -}%
   \fi\fi
  \def\after_setbox{%
    \box_tmp_dim_a\wd\box_tmp
% prepare to revolving:
    \setbox\box_tmp\hbox{%
     \ifx\param_\plus_\kern-\box_tmp_dim_a\fi
     \box\box_tmp
     \ifx\param_\plus_\kern\box_tmp_dim_a\fi}%
% compute dimensions of the box to be revolved:
    \box_tmp_dim_a\ht\box_tmp \advance\box_tmp_dim_a\dp\box_tmp
    \box_tmp_dim_b\ht\box_tmp \box_tmp_dim_c\dp\box_tmp
    \dp\box_tmp0pt \ht\box_tmp\wd\box_tmp \wd\box_tmp\box_tmp_dim_a
% revolve:
    \kern \ifx\param_\plus_ \box_tmp_dim_c \else \box_tmp_dim_b \fi
    \special{ps: gsave currentpoint currentpoint translate
             #190 rotate neg exch neg exch translate}%
    \box\box_tmp
    \special{ps: grestore}%
    \kern -\ifx\param_\plus_ \box_tmp_dim_c \else \box_tmp_dim_b \fi
    \egroup}%
  \afterassignment\jump_setbox\setbox\box_tmp =
}%
\def\revolve{\revolvedir-}
%---------------------------------------------------------------
\def\xflip{% to be followed by an \hbox or a \vbox
  \hbox\bgroup
  \def\after_setbox{%
    \box_tmp_dim_a.5\wd\box_tmp
% prepare to flipping:
   \setbox\box_tmp
     \hbox{\kern-\box_tmp_dim_a \box\box_tmp \kern\box_tmp_dim_a}%
% flip:
   \kern\box_tmp_dim_a
    \special{ps: gsave currentpoint currentpoint translate
         [-1 0 0 1 0 0] concat neg exch neg exch translate}%
    \box\box_tmp
    \special{ps: grestore}%
    \kern-\box_tmp_dim_a
    \egroup}%
  \afterassignment\jump_setbox\setbox\box_tmp =
}%
%---------------------------------------------------------------
\def\yflip{% to be followed by an \hbox or a \vbox
  \hbox\bgroup
  \def\after_setbox{%
    \box_tmp_dim_a\ht\box_tmp \box_tmp_dim_b\dp\box_tmp
    \box_tmp_dim_c\box_tmp_dim_a \advance\box_tmp_dim_c\box_tmp_dim_b
    \box_tmp_dim_c.5\box_tmp_dim_c
% prepare to flipping:
   \setbox\box_tmp\hbox{\vbox{%
     \kern\box_tmp_dim_c\box\box_tmp\kern-\box_tmp_dim_c}}%
% flip:
   \advance\box_tmp_dim_c-\box_tmp_dim_b
   \setbox\box_tmp\hbox{%
     \special{ps: gsave currentpoint currentpoint translate
         [1 0 0 -1 0 0] concat neg exch neg exch translate}%
     \lower\box_tmp_dim_c\box\box_tmp
     \special{ps: grestore}}%
% restore dimensions of the flipped box:
    \ht\box_tmp\box_tmp_dim_a \dp\box_tmp\box_tmp_dim_b
    \box\box_tmp
    \egroup}%
  \afterassignment\jump_setbox\setbox\box_tmp =
}%
%---------------------------------------------------------------
\catcode`\_\undtranscode
%---------------------------------------------------------------
\endinput

bi-cont1.bbl
\bibitem{Bristowe+Crocker}P. \, D. Bristowe and A.\ and G. Crocker,
{\em A computer simulation study of the structures of twin boundaries  in  body-centred cubic crystals},  Philosophical Magazine {\bf 31}, 503(1975)  

\bibitem{copenhagen}D. Rogula,
{\em Non-classical material continua}, in: Theoretical and Applied  Mechanics, IUTAM, (F.I. Niordson, N. Olhoff, eds.),  Elsevier Science Publishers B.V. (North Holland) 1985, p.339  

\bibitem{Rogula+Czerwinska 2001}D. Rogula and T. Lenkowska-Czerwi/nska,
{\em Alternating states in layered high-T\(_c\) superconductors},  Journal of Technical Physics {\bf 42} (2001), p.269  

\bibitem{MS Udine}M. Sztyren,
{\em Solvable nonlocal boundary-value problems}, in: Nonlocal Theory of  Material Media, (D. Rogula, ed.), CISM Courses and Lectures No. 268,  Springer-Verlag 1982  

\bibitem{Kroner 67}E. Kr\"{o}ner,
{\em Elasticity theory of materials with long-range cohesive forces,}  Int. J. Solid Struct. {\bf 3}, 731(1967)  

\bibitem{Walker+Luettmer 96}M. \ and B. Walker and J. Luettmer--Strathmann,
{\em Josephson tunneling in high--$T_c$ superconductors,}   Physical Review B {\bf 54}, 588(1996)  


