# Nonlinear theory of autooscillations of quasiplanar interface during directional solidification

## I INTRODUCTION

A great variety of patterns can form during directional solidification of binary mixture when steady state motion of the planar solid–liquid interface becomes unstable (for a review see, e.g., ). There are several different mechanisms of this instability. One of them, proposed by Mullins and Sekerka  and called the marginal (morphological) instability, has been investigated and discussed intensively for last years . By this mechanism, under appropriate conditions the planar interface loses stability mainly with respect to essentially nonuniform perturbations, which gives rise to cellular interface formation.

Coriell and Sekerka  have shown that departure from local equilibrium at the planar solid/liquid interface due to high solidification rates can cause oscillatory instability of its steady state motion. Merchant and Davis  analyzed this instability in detail taking into account different characteristics of non-equilibrium conditions in a thermodynamically-consistent way. Such an oscillatory instability can manifest itself in self-formation of structures parallel to the moving interface, which have been observed experimentally for Ag–Cu, Al–Cu, and Al–Fe alloys solidified at high rates  (for a review of materials exhibiting formation of this type structures see ). It should be noted that there are also other models predicting oscillatory instability of plane solidification front .

When solidification proceeds by the normal growth mechanism, non-equilibrium effects at the solid/liquid interface become remarkable only at extremely high solidification rates ($``$cm/sec–m/sec) . In this case within the framework of classic approach  departure from local equilibrium at the planar interface due to its rapid motion is characterized by the interface temperature $`T_i(v,C_i)`$ depending on the interface velocity $`v`$ as

$$T_i=T_0+m(v)C_i\frac{v}{\eta }$$
(1)   
and the partition coefficient $`k(v)=C_s/C_i`$ is also treated as a certain function of $`v`$. Hear $`T_0`$ is the melting temperature of the pure material being in thermodynamical equilibrium, $`C_s`$ and $`C_i`$ are the interfacial concentrations of solute in the solid and liquid, respectively, $`\eta `$ is the interface mobility, and $`m(v)`$ is the nonequilibrium slope of the liquidus regarded as a function of the interface velocity $`v`$. Several models for the $`k(v)`$ dependence at high solidification rates have been proposed  and the expression

$$k(v)=\frac{k_e+\beta _0v}{1+\beta _0v}$$
(2)   
is much used in theoretical analysis of pattern formation during rapid directional solidification. Here $`k_e`$ is the equilibrium partition coefficient usually much less than unity, $`k_e1`$, and $`1/\beta _0`$ is the characteristic velocity being typically of the size 5 m/sec . The functions $`m(v)`$ and $`k(v)`$ are not independent but related by the following expression 

$$m=m_e\left\{1+\frac{k_ek[1\mathrm{ln}(k/k_e)]}{1k_e}\right\},$$
(3)   
where $`m_e`$ is the equilibrium value of $`m`$.

In mathematical analysis of directional solidification it is conventionally used the “frozen temperature” approximation. In this approximation the latent heat diffusion is neglected and the thermal properties of the liquid and solid are taken to be equal, which decouples the temperature and solute fields and enables one to represent the temperature distribution $`T(y)`$ in the frame moving at the pulling velocity $`V`$ in the direction $`y`$ in the form

$$T(y)=T_0+Gy,$$
(4)   
where $`G`$ is the temperature gradient at the interface and the origin of the moving coordinate system is chosen so that $`T(y)|_{y=0}=T_0`$. Sivashinsky and Novick-Cohen  have shown that the difference in magnitude of the temperature diffusivity in solid and liquid gives rise to an essentially nonlocal self–interaction of the interface in addition to the ordinary renormalization of the temperature gradient $`G`$ . The effect of latent heat diffusion on rapid solidification has been investigated in detail by Huntley and Davis , and Karma and Sarkissian . It has been found that latent heat diffusion can suppress the longwave oscillatory instability at high solidification rates. However, because of the solute diffusivity $`D`$ in liquid being much less than the temperature diffusivity $`D_\text{T}`$ ( $`D`$ $`D_\text{T}`$) the neighborhood of the zero wave number, where the latent heat effect is considerable, is extremely narrow. The latter enables the temperature $`T_i`$ at the interface to be related to its coordinate $`\zeta (t)`$ along the $`y`$–axis of the moving frame by the expression 

$$T_i(t)=T_0+G\zeta (t)+\frac{L}{c_\text{p}}_{\mathrm{}}^t\frac{dt^{}}{\sqrt{4\pi D_\text{T}(tt^{})}}\frac{d\zeta (t^{})}{dt^{}}.$$
(5)   
Here $`c_\text{p}`$ is the specific heat, $`L`$ is the latent heat of fusion, and $`\zeta (t)`$ is the averaged position of the interface at time $`t`$. When the mean curvature $`K`$ of the interface is not extremely small so that $`Kh_\text{T}1`$ but $`Kh1`$, where $`h_\text{T}D_\text{T}/V`$ and $`hD/V`$ are the characteristic lengths of heat and solute diffusion, respectively, the integral term in expression (5) can be ignorable. In fact, in this case different interface segments of size $`1/K`$ may move independently of one another, thus, their total contribution to the integral term should be small, and the latent heat diffusion controls solely their averaged position $`\zeta (t)`$ along the $`y`$–axis, fixing the value $`\zeta (t)`$. Under these conditions, as follows from (5), the “frozen temperature” approximation holds also for rapid solidification, at least, in semiquantitative analysis. The fact that longitudinal structures forming during rapid solidification are quasi-plane rather than rigorously plane is justified not only by the available experimental data but also numerically.

When crystals grow layer–by–layer, with their surface being practically parallel to one of the singular faces, the solute partitioning can become non-equilibrium at slow rates typically used in technological processes, viz., at $`V10^3`$$`10^2`$ cm/sec . The matter is that the crystal growth at the volumetric solidification rate $`V`$ is attained when atomic steps move along the crystal surface at the mean speed $`V_{\text{step}}V/\theta _{\text{cr}}`$, where $`\theta _{\text{cr}}`$ is the mean angle of interface misorientation from the ideal singular face. For real singular faces forming in crystallization of monocrystals $`\theta _{\text{cr}}10^310^4`$ , thus, $`V_{\text{step}}1100`$ cm/sec for these values of $`V`$. In the layer–by–layer crystallization solute segregation to the solid involves two stages: surface trapping of solute atoms by the moving steps and their following incorporation into the crystal bulk. For such large values of $`V_{\text{step}}`$ the former process becomes essentially non-equilibrium, causing the total partition coefficient $`k(v)`$ to depend on the interface velocity $`v`$. As shown by Voronkov  in this case the $`k(v)`$ dependence can be represented in terms of

$$k(v)=\frac{k_e+k_s\beta _s(v)}{1+\beta _s(v)},$$
(6)   
where the function $`\beta _s(v)`$ is of the form

$$\beta _s(v)=\beta _1v\left[1+\beta _2v+\frac{k_e\beta _2}{k_s\beta _1}\right],$$
(7)   
the constants $`\beta _1`$, $`\beta _2`$ are of order $`10^2`$ sec/cm, $`k_e`$, as before, is the equilibrium partition coefficient, $`k_s`$ is the interface partition coefficient, and typically $`k_ek_s<1`$. The given velocity dependence of the partition coefficient for Al and P in Si has been observed by Voronkov et al. . When the growing interface deviates from the singular face the value $`\theta _{\text{cr}}`$ is ordinary about 1$`^\text{o}`$ or large and the solute partitioning can become non-equilibrium at substantially higher rates of crystallization. A singular crystal surface corresponds to a cusped minimum of the surface energy and, thereby, to an extremely large value of the surface tension . The latter suppresses the marginal instability and practically does not affect the oscillatory instability because at such low crystallization rates the effect of latent heat diffusion should not be strong enough. In this case the non–equilibrium solute partitioning through the induced oscillatory instability can manifest itself in formation of solute bands, where the solute concentration varies only in the growth direction perpendicular to the interface, with the solute bands occurring near singular faces only. Such pattern formation has been found to form during crystallization of InSb with impurities Se and Te in addition to ordinary solute layers caused by temperature fluctuations due to convective instabilities (see  and references therein).

There are a number of works concerned with a nonlinear theory of pattern formation in rapid solidification. Braun and Davis  and Braun et al.  developed a weakly-nonlinear theory and derived the corresponding Ginzburg–Landau equation in the “frozen temperature” approximation. Huntley and Davis  incorporated in this theory latent heat diffusion. Novick-Cohen  obtained a Kuramoto–Sivashinsky equation in the limit $`k1`$. Based on the theory of oscillations with weak damping  Merchant et al.  have studied nonlinear behavior of the zero-wavenumber oscillatory instability in a small neighborhood of the absolute stability point where dissipation processes are not pronounced. The same problem was analyzed numerically by Brattkus and Meiron . Although the found phase paths in the $`v\zeta `$–plane corresponding to the interface oscillations deviate from the ellipsoidal shape substantially their theory describes actually small amplitude interface oscillations as witnessed by the existence of unclosed phase paths in the $`v\zeta `$–plane . The appearance of unclosed phase paths in a small neighborhood of the threshold and the transition from small amplitude (but may be strongly–nonlinear) oscillations to really relaxation oscillations have been analyzed in detail by Baer and Erneux  having considered as an example the well known “French duck solution” of the FitzHugh–Nagumo equation. In  the authors of the present paper actually following the boundary layer model proposed by Langer et al.  for directional solidification have studied strongly nonlinear large amplitude oscillations of planar interface and numerically found anomalous behavior of the corresponding limit cycle. It turns out that although the corresponding nullcline is of the “N–form” the limit cycle goes near only one of its stable branches whereas in the region containing the other stable branch the phase path deviates from it substantially. It should be noted that in the phenomenological model for banded structure formation by Carrad et al.  the phase path has been assumed to follow the nullcline in the conventional way.

In the present paper we develop a theory of strongly nonlinear relaxation oscillations of quasiplanar solid/liquid interface which can occur during directional solidification. The term “quasiplanar” implies that the interface can be treated as planar on spatial scales of order $`hD/V`$ . Since $`D_\text{T}D`$ and so $`h_\text{T}h`$ such an interface does not have to be regarded as plane from the heat transfer standpoint, thus, we may assume that the characteristic curvature radius $`R`$ of this interface meets the inequality $`hRh_\text{T}`$. Besides, as has been mentioned above real longitudinal structures forming during directional solidification are not rigorously plane as well as those observed in numerical experiments . This enables us to make use of the “frozen temperature” approximation, namely, we suppose that the interface temperature $`T_i`$ and the interface position $`\zeta `$ in the coordinate system moving at the pulling velocity $`V`$ are related by expression (4). Dealing with rapid directional solidification this approximation can be justified, at least, in semiquantitive analysis especially if we consider such values of the basic physical parameters that do not correspond to a neighborhood of the nose of the neutral stability curve, where latent heat diffusion affects solidification process most strongly . Besides, in this case to a first approximation in the small parameter $`\zeta /R`$ the nonlocal self–interaction of the interface due to difference in the temperature diffusivities of the solid and liquid vanishes . When crystals grow layer–by–layer approximation (4) is also justified because nonequilibrium partitioning comes into play even at low crystallization rates  at which the latent heat effect is not yet pronounced.

The second assumption adopted in the present analysis is that the interface remains quasi-planar during oscillations. This can be the case when, in particular, the marginal instability is suppressed by large values of surface tension such that the absolute stability limit for the marginal instability lies above that of the oscillatory instability. Such conditions and the corresponding physical parameters have been analyzed in  for rapid solidification. For the layer–by–layer growth the given assumption is justified to a larger degree because a singular crystal surface corresponds to a cusped minimum of the surface energy and, so, to extremely great values of the surface tension . The given assumption allows us to confine ourselves to one–dimensional equations governing solute distribution in the melt.

The third basic assumption to be used is the smallness of the solute partition coefficient $`k(v)`$ for all values of the interface velocity $`v`$ attained during oscillations. For rapid solidification problem this means that we shall consider only the limit of relatively low solidification rates. Concerning with the layer–by–layer growth this assumption does not impose such exacting requirements. Indeed, in this case, as follows from formula (6), $`k(v)1`$ for all values of $`v`$ if $`k_e`$, $`k_s1`$. The latter inequality is justified, e.g., for such crystals as Si and Ge . Besides, in the present analysis we shall ignore the dependence of the liquidus slope $`m`$ on the interface velocity $`v`$ because for small values of the partition coefficient $`k(v)`$ its variations should not lead to appreciable variations of $`m`$ as results from expression (3). The assumption of the partition coefficient $`k(v)`$ being small allows the solute distribution in the melt to be treated as quasistationary and enables us to reduce the original problem to a system of two ordinary differential equations whose solution can be investigated by the standard methods of nonlinear analysis. As will be seen from the results to be obtained below the general properties of interface oscillations are not sensitive to a particular form of the $`k(v)`$ dependence. Keeping the latter in mind, in the following analysis for the sake of definiteness we shall use such a function $`k(v)`$ that generalizes expressions (2) and (6) in a simple way. Besides, for crystals growing layer–by–layer and by the normal growth mechanism the velocity dependence of the interface temperature $`T_i(v)`$ is different. Nevertheless, for the same reason we shall use expression (1).

The paper is organized as follows. In Sec. II we review the equations specifying the directional solidification problem under consideration, briefly discuss linear stability of the steady–state interface motion and recall under what conditions the interface may be treated as quasiplanar. In Sec. III within the framework of the quasistationary approximation the full diffusion problem is reduced to the system of two ordinary differential equations. This section also includes the bifurcation analysis of interface oscillations in order to compare the results obtained in the quasistationary approximation with the previous ones. In Sec. IV we investigate the characteristic properties of the interface relaxation oscillations that occur providing the bifurcation is subcritical and estimate the corresponding parameters of the longitudinal structure forming during solidification.

## II GOVERNING EQUATIONS. THE COEXISTENCE OF DIFFERENT TYPE INSTABILITIES

We confine ourselves to the classic model for directional solidification of a binary mixture . A more introductory and detailed exposition of this process can be found in  and here we only summarize the basic equations governing interface motion and our notation.

The solute concentration $`C`$ in the melt obeys the equation:

$$\frac{C}{t}=D^2C,$$
(8)   
where $`D`$ is the chemical diffusivity in the melt. Far in front of the solid/liquid interface $``$ the concentration $`C`$ tends to $`C_{\mathrm{}}`$

$$CC_{\mathrm{}}\text{as}y\mathrm{}.$$
(9)   
The solute diffusion in the solid is neglected. Local conservation of solute at the interface $``$ leads to the following boundary condition:

$$v_nC_i[1k(v_n)]=D_nC|_{}.$$
(10)   
Here $`v_n`$ is the normal interface velocity, $`C_i`$ is the solute concentration at the interface $``$ on the melt side and $`k(v_n)`$ is the partition coefficient. The velocity dependence of partition coefficient $`k(v_n)`$ due to nonequilibrium solute segregation is specified by the expression

$$k(v)=\frac{k_e+k_s\beta ^{}v}{1+\beta ^{}v},$$
(11)   
where $`\beta ^{}`$ is a given constant. This formula is practically just a simple generalization of expressions (2) and (6).

Keeping in mind the aforesaid in Sec. I we confine our consideration to the “frozen temperature” approximation in which the temperature $`T_i`$ at the interface $``$ and its $`y`$–coordinate $`\zeta `$ in the frame moving at the pulling velocity $`V`$ in the $`y`$–direction are related by the expression

$$T=T_0+G\zeta .$$
(12)   
Here $`T_0`$ is a reference temperature chosen to be equal to the solidification temperature of pure melt ($`C=0`$) and $`G`$ is the imposed temperature gradient.

Following  we suppose that the deviation of the interface temperature $`T_i`$ from the local equilibrium value $`T_{\text{eq}}`$ causes the crystal growth at the velocity:

$$v_n=\eta (T_{\text{eq}}T_i),$$
(13)   
where $`\eta `$ is the interface mobility assumed to be isotropic and constant. The local equilibrium temperature $`T_{eq}`$, solute concentration $`C_i`$ at the interface $``$ and the interface curvature $`H`$ are related by the Gibbs-Thomson condition:

$$T_{\text{eq}}=T_0(1dH)mC_i,$$
(14)   
where $`m>0`$ is the slope of the liquidus line in the binary phase diagram, and $`d`$ is the capillary length which is proportional to the surface tension and accounts for anisotropic effects. The value of $`m`$ may be treated as constant for such $`v_n`$ that $`k(v_n)1`$. In paper  the interface mobility $`\eta `$ has been represented actually as $`\eta =V_0/m`$.

Expressions (8)–(14) form the complete model for the directional solidification process under consideration.

The steady-state solution of equations (8)–(14) describes the motion of the planar interface $``$ at the constant velocity $`v^{\text{st}}=V`$, the uniform distribution of impurity in the solid $`C_{\text{sol}}^{\text{st}}=C_{\mathrm{}}`$, and the stationary solute distribution in the liquid , $`y\zeta `$

$$C^{\text{st}}=C_{\mathrm{}}\left(1+(\frac{1}{k(V)}1)\mathrm{exp}[(y\zeta )\frac{V}{D}]\right).$$
(15)   
Under certain conditions the steady–state motion of the interface loses stability and two different type instabilities, oscillatory and marginal, can occur . The coexistence of these instabilities has been analyzed in general in . However for the purpose of the present paper it will be the best to consider briefly this problem in the limit of small partition coefficient $`(k1)`$.

As follows from the results to be obtained below the amplification rate $`\sigma `$ of a small perturbation $`\xi (x,t)\mathrm{exp}[\frac{V}{D}(\sigma Vt+ipx)]`$ of the planar interface (where $`p`$ is the dimensionless wave number) can be considered to be substantially less than one, $`|\sigma |`$$`1,`$ as long as $`k(v)1`$. In this limit the standard linear stability analysis leads to the following eigenvalue equation:

$$\sigma ^2(1+\mu )+\sigma \mathrm{}(p)+\mathrm{}(p)=0,$$
(16)   
where the functions

$`\mathrm{}(p)`$ $`=`$ $`{\displaystyle \frac{1}{M}}+\mu kVk_V^{}+(2+\gamma +\mu )p^2,`$   
$`\mathrm{}(p)`$ $`=`$ $`{\displaystyle \frac{k}{M}}+[{\displaystyle \frac{1}{M}}+\gamma k1]p^2+\gamma p^4,`$   
and the parameters

$$M=\frac{VmC_{\mathrm{}}}{GDk},\gamma =\frac{T_0dVk}{DmC_{\mathrm{}}},\mu =\frac{Vk}{\eta mC_{\mathrm{}}}.$$

It should be noted that Merchant and Davis  used practically the same parameters where, however, the partition coefficient takes the equilibrium value $`k_e`$.

According to equation (16) the planar interface becomes unstable (i.e. Re$`\sigma >0`$) with respect to the given perturbation when $`\mathrm{}(p)<0`$ or $`\mathrm{}(p)<0`$. The former inequality is associated with the marginal instability, which is caused by the constitutional undercooling and has been thoughtfully analyzed and discussed (for a review see ). The corresponding maximum of the amplification rate $`\sigma _{\text{max}}`$ is attained at $`p=p_{\text{max}}0`$ and, thereby, the neutral curve I, separating stable and unstable (from the marginal instability standpoint) regions in the $`C_{\mathrm{}}V`$–plane, depends on the capillary length $`d`$ (Fig. 1). This instability region is specified by the inequality $`\mathrm{}(p_{\text{max}})<0`$, leading to the expression

$$mC_{\mathrm{}}k\left\{\frac{DG}{V}+\frac{kT_0dV}{D}+\mathrm{\hspace{0.17em}2}\sqrt{T_0Gdk}\right\}.$$
(17)   
The latter inequality is bound up with the oscillatory instability which is described in  and occurs through the dependence of the partition coefficient $`k(v)`$ on the velocity $`v`$. Since the function $`\mathrm{}(p)`$ is increasing the corresponding maximum of amplification rate $`\sigma `$ is attained at $`p=0`$, thus in the “frozen temperature” approximation the oscillatory instability region in the $`C_{\mathrm{}}V`$ \- plane is determined by the inequality

$$mC_{\mathrm{}}\frac{k}{Vk_V^{}}\left\{\frac{DG}{V}+\frac{Vk}{\eta }\right\}$$
(18)   
and is independent of the capillary length $`d`$. The latent heat diffusion can suppress the oscillatory instability for small wave numbers . However due to the heat diffusivity being much larger than the solute one the wave number interval where the latent heat essentially affects perturbation dynamics is narrow enough. Therefore, the condition (18) obtained in the “frozen temperature” model holds true if the latent heat fusion is not extremely great. The equality in (18) specifies the neutral curve II of the oscillatory instability (Fig. 1).

In particular, it follows from (18) that for arbitrary values of the pulling velocity $`V`$ and the imposed temperature gradient $`G`$ all disturbances corresponding to the oscillatory instabiblity decay when $`C_{\mathrm{}}<C^{}`$, where

$$C^{}=\frac{k_e^2}{m\eta (k_sk_e)\beta ^{}}.$$
(19)   
In other words, the inequality $`C_{\mathrm{}}<C^{}`$ is the absolute stability criterion for the oscillatory instability in the limit of small partitition coefficients. It should be noted that this criterion matches the results obtained by Merchant and Davis  in the limit $`V0`$. For $`C_{\mathrm{}}>C^{}`$ the oscillatory instability can occur only in the region $`C_{\mathrm{}}>C_{\text{as}}(\beta )`$ specified by the function $`C_{\text{as}}(\beta )`$ of the dimensionless pulling velocity $`\beta =\beta ^{}V`$

$$C_{\text{as}}(\beta )=C^{}(1+\frac{\beta }{r})^2,$$
(20)   
where the parameter $`r=k_e/k_s`$. In the $`C_{\mathrm{}}V`$–plane the curve $`C_{\mathrm{}}=C_{\text{as}}(\beta )`$ is the boundary of the absolute stability region (line 1 in Fig. 2). In this figure line 2 represents the neutral curve of the oscillatory instability for a fixed value of $`G`$ measured in units of $`G^{}=k_e/(D\eta \beta ^2)`$. In particular, in these terms the neutral curve can be represented as

$$C_{\mathrm{}}=C_{\text{as}}(\beta )\left[1+\frac{G}{G^{}}\frac{1+\beta }{\beta ^2(1+\beta /r)}\right].$$
(21)   
The neutral curve in contrast to the absolute stability curve attains its minimum at $`V=V_{\text{ext}}>0`$. For example, as follows from (21) when $`r1`$ and $`r^2G/G^{}1/r`$ the value $`V_{\text{ext}}`$ belongs to the region $`r^2/\beta ^{}V1/\beta ^{}`$ and so

$$V_{\text{ext}}\frac{1}{\beta ^{}}\left(\frac{rG}{2G^{}}\right).$$
(22)   
As the capillary length increases the neutral curve I for a fixed value of $`G`$ goes upwards in the $`C_{\mathrm{}}V`$–plane and, thus, there can be three characteristic relative positions of curves I and II shown in Fig. 1. When the capillary length is small enough ($`dd_1`$, where $`d_1=d_1(D,G,\eta ,T_0)`$ is such a value of the capillary length that curve I crosses curve II approximately at its minimum point) the region of the oscillatory instability is totally located inside that of the marginal instability. In this case the oscillatory instability of the planar interface seems not to be able to manifest itself in any way because of rapid development of spatially nonuniform interface perturbations (Fig. 1(a)). When $`dd_1`$ there is a certain small region in the $`C_{\mathrm{}}V`$–plane that corresponds solely to the oscillatory instability (Fig. 1(b)). In this case the oscillatory instability can occur. However, as it has been shown numerically , the interface becomes substantially nonplanar during the instability development. In this case the banded structures seem to occur . For $`dd_1`$ the large value of the capillary length, i.e. of the surface tension, suppresses the marginal instability and thereby the interface tends to be planar. It should be noted that in this case (Fig. 1(c)) the developed form of the solidification front will be quasiplanar and solute band formation is expected .

## III QUASISTATIONARY APPROXIMATION

In this section we obtain nonlinear evolution equations for the planar interfacial position $`y=\zeta (t)`$. In general, these equations must contain memory effects; that is a displacement of the interface $``$ from its steady-state position causes a perturbation of the impurity distribution which, in turn, affects the motion of the interface at later times. However, in the limit of small partition coefficient it turns out that the interface velocity varies in time so slowly that it remains practically constant during the time needed for relaxation of the solute distribution in the melt. In other words, it means that all temporal scales of interface dynamics are much larger than the characteristic time ($`\tau _vD/V^2`$) it takes for the steady state solute distribution to form when $`C_i`$ and $`v`$ are fixed.

In the frame of reference moving in the $`y`$–direction at the interface velocity $`v`$ the diffusion equation has the form:

$$\dot{C}=DC_{yy}+vC_y.$$
(23)   
Following  we integrate (23) over $`y`$ from $`\zeta `$ to $`\mathrm{}`$ then taking into account the boundary condition (10) we obtain

$$\frac{d}{dt}_\zeta ^{\mathrm{}}𝑑y[C(y,t)C_{\mathrm{}}]=v[C_{\mathrm{}}k(v)C_i].$$
(24)   
In the quasistationary approximation we can regard the left hand side of equation (23) as a small perturbation and neglect it at lower order of the perturbation technique. In other words, setting the transient term in (23) equal to zero and then solving the obtained equation we find the lower approximation of $`C(y,t)`$: for $`y>\zeta `$

$$C(y,t)=C_{\mathrm{}}+[C_i(t)C_{\mathrm{}}]\mathrm{exp}\left\{\frac{v(t)}{D}(y\zeta )\right\}$$
(25)   
Substituting expression (25) into equation (24) we get

$$D\frac{d}{dt}\left[\frac{C_iC_{\mathrm{}}}{v}\right]=v[C_{\mathrm{}}k(v)C_i].$$
(26)   
Equation (26) along with the equation

$$\frac{dy}{dt}=vV$$
(27)   
and the relationship

$$v=\eta (G\zeta +mC_i)$$
(28)   
resulting from (12)–(14) form a complete description of the quasiplanar interface motion.

The steady state solution of the given system of equations is of the form:

$$v^{\text{st}}=V,C_i^{\text{st}}=\frac{C_{\mathrm{}}}{k(V)},\zeta ^{\text{st}}=\frac{1}{G}\left[\frac{V}{\eta }+\frac{mC_{\mathrm{}}}{k}\right].$$
(29)   
It is convenient for the most of the following analysis to use the dimensionless variables, $`a`$ and $`u`$:

$$a=\frac{\zeta \zeta ^{\text{st}}}{\zeta ^{\text{st}}},u=\frac{vV}{V}$$
(30)   
In terms of this variables, the equations of motion become:

$`\tau {\displaystyle \frac{du}{dt}}`$ $`=`$ $`{\displaystyle \frac{(1+u)^3}{(1+a)}}{\displaystyle \frac{k[u]}{k[0]}}\left(aa_0[u]\right),`$ (31)   
$`\tau {\displaystyle \frac{da}{dt}}`$ $`=`$ $`(\tau \omega )^2u.`$ (32)   
Here we have also introduced the quantities:

$$\tau =\frac{D}{V^2k},(\tau \omega )^2=\frac{1}{kM(1+\mu )},k[u]=k(V_0[1+u]),$$
(33)   
where $`\omega `$ is the frequency of the interface oscillations at the threshold, and the function

$$a_0[u]=\frac{1}{(1+\mu )}\left(\frac{k[0]}{k[u]}+\mu u1\right)+(\tau \omega )^2\frac{k[0]}{k[u]}\frac{u}{(1+u)^2}$$
(34)   
specifies the nullcline of equation (31).

In this way the full diffusion problem of interface dynamics is reduced to the system of ordinary differential equations, which describes the nonlinear oscillations of the quasiplanar interface .

In the present section on the basis of evolution equations (31), (32) we also study the bifurcation mode of the oscillation onset, depending on relation between physical parameters such as $`V,G,\eta `$ etc. In addition, to validate the quasistationary approximation, we compare the form of the interface auto-oscillations obtained in the quasistationary approximation for the supercritical bifurcation and those obtained in the one-dimensional full diffusion problem described in Sec. II.

In order to analyze the bifurcation mode we can make use of the Bogolyubov–Krylov–Mitropolskii technique and the theory of Hopf bifurcation . In this way near the threshold the interface motion is described by the quasiharmonic time dependence of the variable $`a=A(t)\mathrm{cos}(\omega (t)t)`$, where $`\omega (t)\omega `$ and $`A(t)`$ are certain functions of the time $`t`$ being practically constant on the temporal scale $`1/\omega `$. Following the standard procedure  we find that the amplitude $`A(t)`$ of the quasiharmonic oscillations obeys the equations

$$\tau \frac{dA}{dt}=\frac{1}{2}A(ϵ+\alpha A^2),$$
(35)   
where the constants $`ϵ`$ and $`\alpha `$ are specified by the formulae

$$ϵ=a_0^{}[0]=\frac{1}{(1+\mu )k[0]}\left(\mu k[0]+\frac{1}{M}k^{}[0]\right),$$
(36)   
$$\alpha =\frac{1}{8(\tau \omega )^2}\left\{a_0^{\prime \prime \prime }[0]+2a_0^{\prime \prime }[0](3+\frac{k^{}[0]}{k[0]})\right\},$$
(37)   
and the primes on the symbols denote the derivatives of the corresponding functions with respect to the variable $`u`$ taken at the point $`u=0`$. In these terms the instability condition of steady state interface motion takes the form $`a_0^{}[0]<0`$. In other words, for the oscillatory instability to occur the nullcline $`a=a_0[u]`$ should be decreasing in the vicinity of the steady state point ($`u=0`$, $`a=0`$). Taking into account (34) we can rewrite the latter inequality in the form $`\mu +1/(k[0]M)<k^{}[0]/k[0]`$ which as it must is exactly the dimensionless form of inequality (18).

In the following analysis it will be convenient to measure the parameters $`\mu `$ and $`1/(k[0]M)`$ also in units of $`k^{}[0]/k[0]`$. In other words, let us introduce new parameters $`\varphi _g,\varphi _\mu 0`$ defined by the formulae

$$\frac{1}{k[0]M}=\varphi _g\frac{k^{}[0]}{k[0]},\mu =\varphi _\mu \frac{k^{}[0]}{k[0]}.$$
(38)   
In particular, at the points of the neutral curve $`\varphi _g+\varphi _\mu =1`$ and at the absolute stability boundary, $`\mu _{\text{as}}=k^{}[0]/k[0]`$, we get $`\varphi _g=0`$ and $`\varphi _\mu =1`$. Then substituting (34) into (37) and calculating the obtained result at the threshold, $`a_0^{}[0]=0`$, we find an expression for the parameter $`\alpha `$ which for the $`k(v)`$ dependence specified by formula (11) takes the form

$$\alpha =\alpha _0[u]\left(\varphi \varphi _b[u]\right).$$
(39)   
Here the parameter $`\varphi =\varphi _g[0,1]`$ characterizes deviation of the system from the absolute stability boundary, $`\alpha _0[u]`$ is a certain positively definite function of the variable $`u`$, and the function $`\varphi _b[u]`$ is given by the expression

$`\varphi _b[\beta ]={\displaystyle \frac{\beta }{r+\beta }}\left[6{\displaystyle \frac{2\beta }{1+\beta }}{\displaystyle \frac{\beta }{r+\beta }}\right]`$ (41)   
$`\times \left[3+{\displaystyle \frac{(1r)\beta }{(1+\beta )(r+\beta )}}(4{\displaystyle \frac{2\beta }{1+\beta }}{\displaystyle \frac{\beta }{r+\beta }})\right]^1.`$   
According to (41) bifurcation of the oscillatory instability is supercritical if $`\varphi <\varphi _b[\beta ]`$ and subcritical for $`\varphi >\varphi _b[\beta ]`$. As follows from (41), when $`k_ek_s`$, i.e. $`r1`$ there are three limits characterizing different behavior of the function $`\varphi _b[\beta ]`$. In the limit $`\beta r`$ (i.e. $`Vr/\beta ^{}`$) the function $`\varphi _b[\beta ]2\beta /r`$ which actually exactly matching the results obtained by Huntley and Davis  for small pulling velocities in the frozen temperature approximation. When $`r\beta 1`$ (i.e. $`r/\beta ^{}V1/\beta ^{}`$) the value $`\varphi _b[\beta ]5/6`$ and for $`\beta 1`$ ($`V1/\beta ^{}`$) we get that $`\varphi _b[\beta ]1`$ as $`\beta 1`$ and $`\varphi _b[\beta ]>1`$ when $`\beta >1`$. In the latter case, $`\beta >1`$, bifurcation is always supercritical. In the $`C_{\mathrm{}}V`$–plane the curve corresponding to the transition point from supercritical to subcritical bifurcation is specified by the formula $`C_{\text{tr}}=C_{\text{as}}(\beta )/(1\varphi _b[\beta ])`$ (for $`\beta =\beta ^{}V<1`$) and is shown by the curve 3 in Fig. 2. For a fixed value of the temperature gradient $`G`$ the neutral curve in the $`C_{\mathrm{}}V`$–plane crosses the curve $`C_{\mathrm{}}=C_{\text{tr}}(\beta )`$ at a certain point $`(V_c,C_c)`$ and for $`V<V_c`$ as well as $`V>V_c`$ the oscillatory instability onset at the threshold is characterized by subcritical and supercritical bifurcation, respectively. In particular, when $`r^2G/G^{}1/r`$ this point belongs to the region $`r/\beta ^{}V1/\beta ^{}`$ and in this case the values $`V_c`$, $`V_{\text{ext}}`$ are related by the expression $`V_c=(2/5)^{1/3}V_{\text{ext}}`$.

In the region $`V>V_c`$ where bifurcation is supercritical and the coefficient $`\alpha `$ is negative the amplitude of the steady-state quasiharmonic oscillations $`A_0`$ near the threshold is

$$A_0=\left(\frac{ϵ}{\alpha }\right)^{\frac{1}{2}}\left[\frac{k^{}[0]}{k[0]}(\frac{1}{k[0]M}+\mu )\right]^{1/2}.$$
(42)   
For $`V<V_c`$, the coefficient $`\alpha `$ is positive and in order to find the amplitude of the interface oscillations a more complicated analysis is required (this case is considered in Sec. IV).

It should be noted that the results of the present section concerning the bifurcation behavior are in qualitative agreement with those obtained in the previous papers for different physical conditions .

To examine the feasibility of the quasistationary approximation in the supercritical mode we compare the limit cycles of the interface autooscillations as well as the time course $`a(t)`$ and $`u(t)`$ obtained by solving equations (31), (32) and the full one-dimensional problem (8)–(14). For this purpose we have solved numerically the equations mentioned above using the following values of the parameters: $`k_0=0.0125`$,$`k_{\mathrm{}}=0.8`$,$`\beta ^{}V=0.125`$,$`kM=10`$, and $`\mu =0.5`$. In this case the ratio $`k^{}[0]/k[0]0.75`$ and, thus, in the $`C_{\mathrm{}}V`$–plane the equilibrium point of the system is practically located near the boundary of the oscillatory instability region because of $`[k^{}/k(k/M+\mu )]/(k^{}/k)0.2`$. The obtained results are plotted in figures 34 where the dashed and solid lines corresponds to the solutions of equations (31), (32) and the full diffusion model, respectively.

As seen from Fig. 3 the quasistationary approximation gives an adequate description of the interface oscillations at least in the supercritical mode. Fig. 4 also demonstrates that even in the vicinity of the instability region boundary the nonlinear effects play a significant role in the interface oscillations and the limit cycle substantially deviates from the elliptic shape and the time dependence $`u(t)`$, $`a(t)`$ are practically of the spikewise form. Nevertheless, the period $`T`$ of these interface autooscillations can be estimated by the expression $`T=2\pi \omega ^1=2\pi \tau [(1+\mu )kM]^{1/2}`$ which is rigorously justified for quasiharmonic oscillations only. Fig. 5 shows the resulting impurity distribution in the solid $`C_{\text{sol}}(y)`$. According to Fig. 5 and as would be expected the spatial period of the growing superlattice is of order

$$HTV2\pi \left[(1+\mu )\frac{DmC_{\mathrm{}}}{Vk^2G}\right]^{1/2}.$$
(43)   
## IV RELAXATION OSCILLATIONS OF INTERFACE

As shown in Sec. III the subcritical bifurcation corresponds to small values of the pulling velocity, namely, $`V<1/\beta ^{}`$ ($`\beta <1`$). Besides, the asymptotic behavior of the curve separating in the $`C_{\mathrm{}}V`$–plane the regions of the subcritical and supercritical bifurcation is different for $`\beta r`$ and $`r\beta 1`$ (when $`r=k_e/k_s1`$). So for $`\beta <1`$ the interface oscillations may be expected to be relaxation, with their main features depending substantially on the ratio of $`\beta `$ and $`r`$. In agreement with the results to be obtained below there are two cases of strongly–nonlinear dynamics of the interface motion: weakly– and strongly–dissipative oscillations. Weakly–dissipative oscillations (however of large amplitude with respect to $`u`$) can occur when

$$\beta r\text{and}\frac{\beta }{r}(1\varphi _\mu ),\varphi _g\left(\frac{\beta }{r}\right)^{1/3}.$$
(44)   
whereas the conditions

$$\beta r\text{and}\left(\frac{\beta }{r}\right)^{1/3}(1\varphi _\mu ),\varphi _g$$
(46)   
or

$$r\beta 1(\text{for }r=k_e/k_s1)$$
(47)   
will result in strongly–dissipative oscillations. It should be noted that in the strict sense only strongly–dissipative oscillations are relaxation, although in the two cases the velocity amplitude $`u`$ can attain large values much greater than unity. In this section we analyze strongly–nonlinear dynamics of the interface oscillations in the given cases individually.

Weakly–dissipative oscillations of crystal interface that can occur during directional solidification has been analyzed in detail by Merchant et al.  and Brattkus and Meiron  within the framework of the classic model for rapid solidification similar to one used in the present analysis. So we only briefly consider the first limit (44) in order to complete the analysis of strongly–nonlinear dynamics of the given model.

In this case the ratio

$$\frac{k^{}[0]}{k[0]}=\frac{(1r)\beta }{(1+\beta )(r+\beta )}\frac{\beta }{r}1$$

is a small parameter, $`\mu 1`$, and the value

$$(\tau \omega )^2\varphi _g\frac{\beta }{r}.$$

According to Sec. III, for $`\beta r`$ the bifurcation is subcritical when $`\varphi _g>2\beta /r`$ and at the threshold $`\varphi _g+\varphi _\mu =1`$. So in order to analyze relaxation oscillations of the interface in the given limit we may confine ourselves to the case $`1\varphi _\mu \varphi _g`$ and $`\varphi _g\beta /r`$. We note that the two relations have actually led us to the left–hand side of the latter conditions of limit (44). Its right–hand side enables us to treat the ratio $`a_0[u]/a`$ in equation (31) as a small parameter for the interface oscillations that will occur under these conditions. The matter is that the characteristic value of the amplitude $`\overline{a}`$ of the variable $`a`$ turns out to be about $`\overline{a}\tau \omega `$, so $`\overline{a}(\varphi _g\beta /r)^{1/2}`$ and the maximum $`u_{\text{max}}`$ of the interface velocity $`u`$ attained during oscillations can be estimated as $`u_{\text{max}}\stackrel{~}{u}_{\text{min}}`$, where $`\stackrel{~}{u}_{\text{min}}`$ is the $`u`$–coordinate of the point $`(\stackrel{~}{u}_{\text{min}},\stackrel{~}{a}_{\text{min}})`$ at which the curve $`a_0[u]`$ attains its extremum (minimum) in the region $`u>0`$ (Fig. 6). In this case expression (34) may be rewritten in the form

$$a_0[u]\frac{\beta }{r}u\left[(1\varphi _\mu )+\frac{\beta }{r}u+\varphi _g\frac{1}{(1+u)^2}\right],$$
(48)   
whence for $`\varphi _g1\varphi _\mu `$ we get that

$$\stackrel{~}{u}_{\text{min}}\frac{r}{2\beta }(1\varphi _\mu )\frac{r}{\beta }\varphi _g1,$$
(50)   
$$\stackrel{~}{a}_{\text{min}}\frac{1}{4}(1\varphi _\mu )^2\varphi _g^2,$$
(51)   
so, by virtue of (44), we obtain that $`|a_0[u]||\stackrel{~}{a}_{\text{min}}|\overline{a}`$.

At lower order in the small parameter $`a_0[u]/a`$ the system of equations (31), (32) is conservative with the “energy” (the first integral)

$$(u,a)=\frac{(\tau \omega )^2}{2}\frac{u^2}{(1+u)^2}+\frac{1}{2}a^2$$
(52)   
When deriving expression (52) we have also taken into account that for the given values of the parameters $`k[u]k[0]`$ and $`a1`$. The remaining term in equation (31) proportional to $`a[u]`$ causes time variation of the “energy”, namely,

$$\tau \frac{d}{dt}=(\tau \omega )^2a_0[u]u.$$
(53)   
In the case under consideration the system dynamics may be described in terms of fast motion along a phase path, specified by the equation $`(u,a)=h`$ for a fixed value of $`h`$, and slow time variations in the “energy” $`(u,a)`$. The geometry of the phase paths in the $`ua`$–plane is shown in Fig. 6. As seen in Fig. 6 there are two types of the phase paths, closed and unclosed, which correspond to $`h<\frac{1}{2}(\tau \omega )^2`$ and $`h>\frac{1}{2}(\tau \omega )^2`$, respectively, and under the adopted assumptions solely the closed paths are meaningful.

Steady–state oscillations matches such a value $`h_{\text{st}}`$ for which the right–hand side of equation (53) averaged over a single period of the system motion is equal to zero. Since the right–hand side of equation (53) depends on the variable $`u`$ only, at lower approximation we may average it assuming the system to move strictly along the phase for a certain fixed value of $`h`$. In this way from (53) we get

$$\tau \frac{dh}{dt}=\frac{2\tau }{T(h)}\underset{u_{\text{min}}(h)}{\overset{u_{\text{max}}(h)}{}}𝑑u\frac{a_+(h,u)}{u}a_0[u],$$
(54)   
where $`T(h)`$ is the period of the system motion along the phase path corresponding to $`=h`$, $`u_{\text{min }}(h)`$ and $`u_{\text{max }}(h)`$ are the minimum and maximum of the variable $`u`$ attained during the motion along this phase path, and $`a_+(h,u)`$ is the positive solution of the equation $`(u,a)=h`$, namely,

$$u_{\text{min}}(h)=\frac{(2h)^{1/2}}{\tau \omega +(2h)^{1/2}},u_{\text{max}}(h)=\frac{(2h)^{1/2}}{\tau \omega (2h)^{1/2}}$$
(55)   
and

$$a_+(h,u)=(\tau \omega )\left[\frac{2h}{(\tau \omega )^2}\left(\frac{u}{1+u}\right)^2\right]^{1/2}.$$
(56)   
Substituting (55) and (56) into (54) and integrating over $`u`$ we obtain

$`\tau {\displaystyle \frac{dh}{dt}}=\pi (\tau \omega )^2{\displaystyle \frac{\beta \tau }{rT(h)}}{\displaystyle \frac{2h}{(\tau \omega )^2}}`$ (58)   
$`\times \left\{(1\varphi _\mu )\mathrm{\Lambda }_1\left[{\displaystyle \frac{2h}{(\tau \omega )^2}}\right]{\displaystyle \frac{\beta }{r}}\mathrm{\Lambda }_2\left[{\displaystyle \frac{2h}{(\tau \omega )^2}}\right]\varphi _g\right\},`$   
where the functions

$$\mathrm{\Lambda }_1[x]=\frac{2}{\sqrt{1x}(1+\sqrt{1x})},$$
(59)   
$$\mathrm{\Lambda }_2[x]=\frac{2x(1+2\sqrt{1x})}{(1x)^{3/2}(1+\sqrt{1x})^2}.$$
(60)   
In agreement with the results obtained in Sec. III from (58) we can see that the oscillatory instability occurs when $`\varphi _\mu +\varphi _g<1`$ and is subcritical for $`\varphi _\mu +2\beta /r<1`$. As follows from (58) the dependence of the “energy” $`h_{\text{st}}`$ of steady state oscillations on the physical parameters is specified by the expression

$$\varphi _g=(1\varphi _\mu )\mathrm{\Lambda }_1\left[\frac{2h_{\text{st}}}{(\tau \omega )^2}\right]\frac{\beta }{r}\mathrm{\Lambda }_2\left[\frac{2h_{\text{st}}}{(\tau \omega )^2}\right].$$
(61)   
In particular, for $`\beta /r1\varphi _\mu (\beta /r)^{3/2}`$ and $`\varphi _g[(1\varphi _\mu )(r/\beta )^{1/3}]^{3/2}`$ the value

$$h_{\text{st}}\frac{(\tau \omega )^2}{2}\left[1\frac{\beta }{r(1\varphi _\mu )}\right].$$
(62)   
Besides, according to (52) the maximum $`a_{\text{max}}`$ and the minimum $`a_{\text{min}}`$ attained by the variable $`a`$ during oscillations are $`a_{\text{max}}a_{\text{max}}(2h_{\text{st}})^{1/2}`$. Whence for the given values of $`\varphi _g`$ and $`\varphi _\mu `$ we obtain

$$\begin{array}{cc}u_{\text{min}}\frac{1}{2},\hfill & a_{\text{min}}\tau \omega ,\hfill \\ u_{\text{max}}\frac{2r}{\beta }(1\varphi _\mu ),\hfill & a_{\text{max}}\tau \omega .\hfill \end{array}$$
(63)   
The period $`T`$ of these oscillations is practically determined by the time $`\tau _s^{}`$ it takes for the system to pass through the region $`u<0`$. Equation (32) enables us to estimate this time as $`\tau _s^{}\tau a_{\text{max}}/(\tau \omega )^2`$ whence we get

$$T\frac{\tau }{\tau \omega }.$$
(64)   
In the case under consideration the characteristic relative position of the nullcline $`a_0[u]`$ and the limit cycle of the oscillations are demonstrated in Fig. 6.

In limits (46) and (47) the interface oscillations are strongly–dissipative. The corresponding phase paths as a whole cannot be described by level curves of any function similar to the energy $`(u,a)`$. In particular, in the region $`u1`$ these phase paths may be roughly related to the unclosed level lines of the energy $`(u,a)`$. In order to analyze the interface oscillations in this case we may make use of the classic theory of relaxation oscillations . First, we consider in detail limit (46). Under these conditions, as it can be shown directly from the system of equations (31), (32), the phase path going through a point that is not located in the region $`|u+1|1`$ or in a small neighborhood of the nullcline $`a=a_0[u]`$ forms practically a horizontal line, at least, in the vicinity of this point. In other words, at such a point the value $`du/dt`$ is large in comparison with $`da/dt`$ and the quantity $`\tau \omega `$ may be treated as a small parameter. Therefore, the interface dynamics should include fast motion governed by equation (31) for fixed $`a`$ until the phase path reaches a small neighborhood of the nullcline $`a=a_0[u]`$, where the left–hand side of equation (31) tends to zero, or a small neighborhood of the boundary $`u=1`$, which formally can be also regarded as the other branch of nullcline of equation (31). In the two latter regions the interface motion is slow.

The nullcline $`a=a_0[u]`$ (shown by curve 1 in Fig. 7) looks like an “N” and the stationary point $`(u=0,a=0)`$ is unstable when it belongs to the decreasing segment of the nullcline. So according to the classic theory of relaxation oscillations  one could expect that the limit cycle will be of the form presented by the dashed line in Fig. 7. Indeed, based on the system of equations (31), (32) we can show that the increasing segments I, III of the nullcline $`a=a_0[u]`$ are formally attractive for the phase path whereas the increasing one is (segment II) is repulsive. The main characteristics of such a limit cycle are directly determined by the form of the nullcline. In particular, the minimum and maximum attained by the variable $`a`$ during oscillations are approximately equal to $`\stackrel{~}{a}_{\text{min}}`$ and $`\stackrel{~}{a}_{\text{max}}`$, i.e. the $`a`$–coordinates of the extremum points ($`\stackrel{~}{u}_{\text{min}},\stackrel{~}{a}_{\text{min}}`$), ($`\stackrel{~}{u}_{\text{max}},\stackrel{~}{a}_{\text{max}}`$) of the nullcline. In the case under consideration the curve $`a_0[u]`$ is specified by expression (48) and the last term on its right–hand side is ignorable for $`u1`$. Whence we find that the quantities $`\stackrel{~}{a}_{\text{min}},\stackrel{~}{u}_{\text{min}}`$ can be estimated by expressions (63) and for $`\stackrel{~}{a}_{\text{max}},\stackrel{~}{u}_{\text{max}}`$ we obtain

$$|\stackrel{~}{u}_{\text{max}}|10<\stackrel{~}{a}_{\text{max}}\frac{(1\varphi _\mu )\beta }{r}$$
(65)   
However, the right–hand side of equation (31) possesses a certain peculiarity in the region $`|u+1|1`$, because it tends to zero as $`u1`$. Therefore, if the phase path goes into this region the time scale hierarchy will be changed and, thus, the phase path will not be able to follow the segment of the nullcline $`a=a_0[u]`$ going through this region, and the maximum $`a_{\text{max}}`$ attained by the variable $`a`$ during oscillations will be not determined by the extremum point ($`\stackrel{~}{u}_{\text{max}},\stackrel{~}{a}_{\text{max}}`$). This is the case in limit (46) because for such values of the physical parameters the nullcline $`a=a_0[u]`$ reaches a small neighborhood of the line $`u=1`$ in the region $`a>0`$ and the phase path has to go into the region $`|u+1|1`$ returning from large values of $`u`$. Therefore, under the given conditions the limit cycle of the interface oscillations will involve the conventional segments going along the horizontal lines $`a=a_{\text{max}}`$, $`a=\stackrel{~}{a}_{\text{min}}`$, along segment III of the nullcline $`a=a_0[u]`$ from $`a_{\text{max}}`$ to $`\stackrel{~}{a}_{\text{min}}`$ in addition to a certain anomalous segment located in the region $`|u+1|1`$ which joins the two horizontal segments.

In order to complete the limit cycle construction we need to analyze the system motion in the region $`|u+1|1`$ by solving directly the system of equations (31), (32). Dividing equation (31) by equation (32), setting $`u=1`$ except for the terms containing the cofactor $`(1+u)`$, and taking into account that in the given case $`a1`$ we obtain the following equation describing the phase path in the region $`|1+u|1`$

$$\frac{du}{da}=(1+u)^3\left[\frac{r}{\beta \varphi _g}a+\frac{1}{(1+u)^2}\frac{1\varphi _\mu }{\varphi _g}\right]$$
(66)   
subject to the formal “initial” condition

$$u\mathrm{}\text{as}a\stackrel{~}{a}_{\text{min}}+0.$$
(67)   
Condition (67) reflects the fact that the limit cycle segment located in the region $`|1+u|1`$ originates from the phase path going into this region practically along the horizontal line $`a=\stackrel{~}{a}_{\text{min}}`$. To second order in $`a`$ the solution of equation (66) meeting (67) is of the form

$$\frac{1}{(1+u)^2}\frac{r}{\beta \varphi _g}(a\stackrel{~}{a}_{\text{min}})(a\stackrel{~}{a}_{\text{min}}).$$
(68)   
As it must $`u\mathrm{}`$ as $`a\stackrel{~}{a}_{\text{min}}+0`$, the second point where the value $`u`$ tends to infinity is $`a=\stackrel{~}{a}_{\text{min}}`$. Besides, according to (68), values of $`u`$ such that $`(1+u)1`$ correspond to values of $`a`$ belonging to small neighborhoods of the points $`\stackrel{~}{a}_{\text{min}}`$ and $`\stackrel{~}{a}_{\text{min}}`$, i.e. $`|a+\stackrel{~}{a}_{\text{min}}||\stackrel{~}{a}_{\text{min}}|`$ or $`|a\stackrel{~}{a}_{\text{min}}||\stackrel{~}{a}_{\text{min}}|`$ because in the given case the ratio $`r\stackrel{~}{a}_{\text{min}}^2/(\beta \varphi _g)`$ is a large parameter. Whence it follows that the desired value of $`a_{\text{max}}`$ can be estimated as $`a_{\text{max}}\stackrel{~}{a}_{\text{min}}`$ i.e.

$$a_{\text{max}}\frac{(1\varphi _\mu )^2}{4}.$$
(69)   
It should be noted that $`a_{\text{max}}\stackrel{~}{a}_{\text{max}}`$, so the limit cycle goes over the extremum point ($`\stackrel{~}{u}_{\text{max}},\stackrel{~}{a}_{\text{max}}`$) at a remarkable distance. In addition, formula (68) shows that the minimum $`u_{\text{min}}`$ attained by the variable $`u`$ during oscillations is about $`u_{\text{min}}1+[(\beta \varphi _g)/(r\stackrel{~}{a}_{\text{min}}^2)]^{1/2}`$, thus

$$u_{\text{min}}1+\left[\frac{4\beta \varphi _g}{r(1\varphi _\mu )^4}\right]^{1/2}.$$
(70)   
The maximum $`u_{\text{max}}`$ attained during oscillations is determined by the intersection point of the horizontal line $`a=a_{\text{max}}`$ and the nullcline $`a=a_0[u]`$, whence we get

$$u_{\text{max}}\frac{(1+\sqrt{2})(1\varphi _\mu )r}{2\beta }.$$
(71)   
Let us now estimate the characteristic time scales of the interface oscillations. According to the classic theory of relaxation oscillations  the fast motion along the horizontal lines $`a=\stackrel{~}{a}_{\text{min}}`$ and $`a=a_{\text{max}}`$ governed by equation (31) is characterized by the time scale $`\tau _f\tau /a_{\text{max}}`$. The time during which the system moves slowly inside the region $`|1+u|1`$ practically in the $`a`$–direction from the point $`\stackrel{~}{a}_{\text{min}}`$ to $`a_{\text{max}}`$ is about $`\tau _s^{}2a_{\text{max}}\tau /(\tau \omega )^2\tau (1\varphi _\mu )^2r/(\beta \varphi _g)\tau _f`$. The motion along segment III of the nullcline $`a=a_0[u]`$ is actually governed by equation (32) where we may roughly set $`uu_{\text{max}}`$. In this way we find that for this stage of the interface motion the characteristic time can be estimated as $`\tau _s^+\tau _s^{}/u_{\text{max}}\tau (1\varphi _\mu )/\varphi _g\tau `$. It should be pointed out that in spite of $`\tau _s^+\tau _f\tau /a_{\text{max}}`$ the motion along the nullcline $`a_0[u]`$ in the region $`uu_{\text{max}}`$ is slow in comparison with the fast motion along the horizontal lines $`a=a_{\text{max}}`$ or $`a=a_{\text{min}}`$ because in this region the time scale of the fast motion is $`\tau _f/u_{\text{max}}^3`$ rather than $`\tau _f`$. In the given case the period $`T`$ of the system motion is determined, as before, by the value $`\tau _s^{}`$, thus

$$T\tau \frac{r}{\beta }\frac{(1\varphi _\tau )^2}{\varphi _g}$$
(72)   
Let us now consider the interface oscillations in limit (47). In addition we shall assume that $`\varphi _\mu \varphi _{\mu b}=1/6`$ (Sec. III), because for $`\varphi _\mu <\varphi _{\mu b}`$ the oscillation onset is subcritical. In this case $`\mu 1`$, the $`k[u]`$–dependence and the nullcline $`a_0[u]`$ may be approximately specified by the expressions

$$k[u]k[0](1+u)$$
(73)   
and

$$a_0[u]\frac{u}{1+u}+\varphi _g\frac{u}{(1+u)^3}+\varphi _\mu u.$$
(74)   
For $`\varphi _g1`$ the system of equations (31), (32) does not contain actually a small parameter. Nevertheless, the minimum $`\stackrel{~}{a}_{\text{min}}`$ attained by the nullcline $`a_0[u]`$ in the region $`u>0`$ is located in the immediate vicinity of the boundary $`a=1`$, i.e. $`\stackrel{~}{a}_{\text{min}}`$ $`1`$ and matches the value $`\stackrel{~}{u}_{\text{min}}1/(\varphi _\mu )^{1/2}1`$. Besides, the increasing segment of the nullcline $`a_0[u]`$ (for $`u>0`$) corresponds to large values of the variable $`u`$. Therefore, in the given case the variable $`u`$ attains large values during interface oscillations and the limit cycle comes close to the line $`a=1`$. Since the right–hand side of equation (31) depends on $`u`$ and $`a`$ as $`(1+u)^4/(1+a)`$ this behavior of the limit cycle causes the interface oscillations to be relaxation and the limit cycle contains segments of the fast motion that approximately are parallel to the $`u`$–axis, a segment following the increasing part of the nullcline $`a=a_0[u]`$ in the region $`u1`$, and a segment located in the region $`|u+1|1`$. The latter one, however, also approximately follows the corresponding increasing part of the nullcline $`a=a_0[u]`$ as results from numerical analysis although the bifurcation is subcritical for the given physical parameters. This characteristics explains the fact that the region of the parameter $`\varphi _g`$ where the interface motion is linearly–stable and the interface oscillations with a finite amplitude can occur is sufficiently narrow. In particular, as formally $`\varphi _\mu 0`$ this region is determined by the inequality $`1<\varphi _g<\varphi _{gc}1.05`$.

Inside the instability region far from the neutral curve, i.e. for $`\varphi _g1`$, the limit cycle again will reach the region $`|u+1|1`$, thus the corresponding segment of the limit cycle will not follows the increasing part of the nullcline $`a=a_0[u]`$ and take a form similar to that matching limit (46). This conclusion has been obtained by solving numerically the system of equation (31), (32) in the region $`|u+1|1`$ for $`\varphi _g1`$ where it can be rewritten in the parameterless form

$`{\displaystyle \frac{u_p}{t}}`$ $`=`$ $`{\displaystyle \frac{u_p^4}{a}}_p\left[a_p{\displaystyle \frac{1}{u}}_p(1{\displaystyle \frac{1}{u_p^2}})\right],`$ (76)   
$`{\displaystyle \frac{a_p}{t}}`$ $`=`$ $`1,`$ (77)   
where $`u_p=(u+1)/\varphi _g^{1/2}`$, $`a_p=(a+1)\varphi _g^{1/2}`$, and the time $`t`$ is measured in units of $`\tau /\varphi _g^{3/2}`$. In the given region the phase path substantially deviates from the nullcline and goes over its extremum point ($`\stackrel{~}{u}_{\text{max}},\stackrel{~}{a}_{\text{max}}`$) at a certain distance. So also in the case under consideration the limit cycle of the interface oscillations and one formally predicted by the classic theory of relaxation oscillations  differ significantly in form. The period of such oscillations can be roughly estimated as

$$T\frac{\tau }{\varphi _g^{3/2}}.$$
(78)   
Concluding the present section we note that, first, in all the cases considered hear the interface oscillations may be regarded as spikewise. The matter is that the interface dynamics described in terms of the variables $`a`$ and $`u`$, i.e. by the system of equations (31), (32) is characterized by self–acceleration in the region $`u1`$, which is reflected in the existence of the cofactor $`(1+u)^3`$ in equation (31). So the interface dynamics is characterized, at least, by two time scales, $`\tau ^{},\tau ^+`$, corresponding to the motion in the regions $`|u|1`$ and $`u1`$, respectively, with $`\tau ^{}\tau ^+`$. Therefore, the time course of these oscillations contains spikes of duration $`\tau ^+`$ formed during motion in the region $`u1`$. Second, from the standpoint of the crystal growth theory it is important to analyze the main characteristics of impurity distribution in the growing crystal. So we also estimate the spatial period $`H`$ of the impurity distribution $`C_{\text{sol}}(y)`$ occurring during strongly–dissipative interface oscillations. By definition

$$T_y=V_0^T𝑑t(1+u).$$
(79)   
For example, in limit (46) the main contribution to the value of integral (79) is due to the system motion along segment III of the nullcline. Then substituting the corresponding estimates into (79) we get $`HVu_{\text{max}}\tau _s^+`$, so

$$H\frac{D}{Vk(V)}\frac{r(1\varphi _\mu )^2}{\beta \varphi _g}$$
(80)   
In case (47) for $`\varphi _g1`$ and $`\varphi _\mu 1`$ the spatial period is about

$$H\frac{D}{Vk(V)}$$
(81)   
The comparison of (80) and (81) shows that expression (81) can be treated as a rough estimate of the spatial period of the impurity distribution occurring during strongly–dissipative interface oscillations. Keeping in mind the layer–by–layer crystal growth we set $`D510^5`$ cm<sup>2</sup>/sec, $`V10^2`$ cm/sec and $`k(V)0.1`$ for such values of $`V`$. Then from (81) we get $`H500`$ $`\mu `$m, which coincides in order with the characteristic length of impurity microinhomogeneities observed in real crystals .

In addition, we consider some characteristic features of the $`C_{\text{sol}}(y)`$ dependence. Expressions (28), (30) lead us to the following relationship

$`{\displaystyle \frac{C_{\text{sol}}}{C_{\mathrm{}}}}`$ $`=`$ $`1+{\displaystyle \frac{k[u]}{k[0]}}(\mu +1)(aa_0[u])+(\tau \omega )^2{\displaystyle \frac{u}{(1+u)^2}}`$ (83)   
$`={\displaystyle \frac{k[u]}{k[0]}}[1+(1+\mu )a\mu u].`$   
During the system motion practically along segment III of the nullcline $`a_0[u]`$ the concentration $`C_{\text{sol}}`$ is actually equal to $`C_{\mathrm{}}`$ because the right hand side of (83) differs from one by the value $`(\tau \omega )^2u/(1+u)^21`$ for $`u1`$. As seen from expression (83) in cases where $`\mu 1`$ the impurity concentration in the solid near the interface attains its minimum $`C_{\text{sol}}^{\text{min}}`$ and maximum $`C_{\text{sol}}^{\text{max}}`$ when in the $`ua`$–plane the system reaches the points with the coordinates $`aa_{\text{min}},uu_{\text{min}}`$ and $`aa_{\text{max}},uu_{\text{max}}`$, respectively. So the minimum and the maximum of the impurity distribution $`C_{\text{sol}}(y)`$ can be estimated as

$`C_{\text{sol}}^{\text{min}}`$ $``$ $`C_{\mathrm{}}{\displaystyle \frac{k[u_{\text{min}}]}{k[0]}}(1+a_{\text{min}}),`$   
$`C_{\text{sol}}^{\text{max}}`$ $``$ $`C_{\mathrm{}}{\displaystyle \frac{k[u_{\text{max}}]}{k[0]}}(1+a_{\text{max}}).`$   
In order to analyze the feasibility of the quasistationary approximation for the interface relaxation oscillations we have also solved the system of equations (31)-(32) as well as the full one-dimensional model (8)-(14) numerically. We used the following values of the parameters: $`k_0=0.0125,k_{\mathrm{}}=0.8,\beta ^{}V=0.125,Mk=10`$, and $`\mu =0.2`$, which correspond to the subcritical bifurcation and are not too far from the critical point. The obtained results are presented in Fig. 8–Fig. 10.

It should be expected that the quasistationary approximation may be violated when the interface velocity becomes small enough: $`(1+u)1`$. Nevertheless, as seen from these figures, the shapes of the limit cycle constructed by different methods qualitatively coincides with each other. Therefore, first, based on the quasistationary approximation one can describe the main characteristics of oscillatory zoning. Second, the estimates for the amplitude and spacing of the growing superlattice can be used as a first approximation for the real processes.

## ACKNOWLEDGMENTS

The authors thank Professor V.V. Voronkov for interest and useful discussions. The research described in the present paper was made possible in part by Grants No. U1I000, U1I200 from the International Science Foundation.