Localization in adiabatic shear flow
via geometric theory of singular perturbations
Abstract
We study localization occurring during high speed shear deformations of metals leading to the formation of shear bands. The localization instability results from the competition between Hadamard instability (caused by softening response) and the stabilizing effects of strainrate hardening. We consider a hyperbolicparabolic system that expresses the above mechanism and construct selfsimilar solutions of localizing type that arise as the outcome of the above competition. The existence of selfsimilar solutions is turned, via a series of transformations, into a problem of constructing a heteroclinic orbit for an induced dynamical system. The dynamical system is four dimensional but has a fastslow structure with respect to a small parameter capturing the strength of strainrate hardening. Geometric singular perturbation theory is applied to construct the heteroclinic orbit as a transversal intersection of two invariant manifolds in the phase space.
1 Introduction
Shear bands are narrow zones of intensely localized shear that are formed during the high speed plastic deformations of metals [31, 2, 30]. They often precede rupture and are one of the striking instances of material instability leading to failure. Considerable attention has been devoted to the problem of shear band formation in both the mechanics and the applied mathematics literature, and section 2 is devoted to a presentation of the problem and a quick derivation of the hyperbolicparabolic system
\displaystyle v_{t}  \displaystyle=\big{(}\theta^{\alpha}\gamma^{m}v_{x}^{n}\big{)}_{x},  (1)  
\displaystyle\gamma_{t}  \displaystyle=v_{x},  
\displaystyle\theta_{t}  \displaystyle=\theta^{\alpha}\gamma^{m}(v_{x})^{n+1}. 
The system describes the plastic shearing deformation of a specimen based on conservation of momentum and energy using a model in thermoviscoplasticity:
\sigma=\theta^{\alpha}\gamma^{m}u^{n}\,,\qquad\mbox{ where \quad$u:=\gamma_{t% }=v_{x}$ }  (2) 
Equation (2) is viewed as a yield stress or a plastic flow rule, with the parameters \alpha, m and n>0 describing respectively the degree of thermal softening, strain hardening and strainrate sensistivity. We refer to section 2 for a derivation of (1) and a review of earlier work useful in understanding the localization problem and its relevance to the present study; references [2, 24, 30, 17] can be consulted for further information on the mechanical aspects of the model.
The model (1) admits a special class of timedependent solutions describing uniform shear (see (17)) and the problem of shear band formation is initially posed as a problem of stability for the uniform shearing solutions. As these are timedependent, it leads to analysis of nonautonomous problems and presents challenges even for linearized stability. We refer to section 2 for information on this aspect of the problem. Here, we focus on the regime of linearized instability and pose the problem of understanding the behavior in the nonlinear regime. A conjecture on the threshold of instability is offered by the asymptotic analysis in [17], devising an effective equation that changes type along a threshold from forward to backward parabolic. It leads one to expect instability when \alpha+m+n changes sign from positive to negative value.
It is expedient to reformulate the problem (1), in terms of the variables (u,\gamma,\theta), as a parabolic system where the diffusion coefficient is controlled by ordinary differential equations
\displaystyle u_{t}  \displaystyle=\big{(}\theta^{\alpha}\gamma^{m}u^{n}\big{)}_{xx},  (3)  
\displaystyle\gamma_{t}  \displaystyle=u,  
\displaystyle\theta_{t}  \displaystyle=\theta^{\alpha}\gamma^{m}u^{n+1}\,. 
The systems (1) or (3) are considered for x\in\mathbb{R}, t>0. The goal of this work is to construct a class of selfsimilar solutions for systems (1) (or (3)) of the form
\displaystyle\gamma(t,x)  \displaystyle=t^{a}\Gamma\big{(}x\,t^{\lambda}\big{)},  \displaystyle v(t,x)  \displaystyle=t^{b}V\big{(}x\,t^{\lambda}\big{)},  \displaystyle\theta(t,x)  \displaystyle=t^{c}\Theta\big{(}x\,t^{\lambda}\big{)},  (4)  
\displaystyle{\sigma}(t,x)  \displaystyle=t^{d}\Sigma\big{(}x\,t^{\lambda}\big{)},  \displaystyle u(t,x)  \displaystyle=t^{a1}U\big{(}x\,t^{\lambda}\big{)}, 
where \lambda>0 and the parameters (\alpha,m,n) take values in the expected instability regime \alpha+m+n<0. Usually parabolic systems (such as (3)) admit diffusing self similar solutions constant on lines \xi=\frac{x}{t^{\rho}}. By insisting on \lambda>0, the solutions (4) will propagate information on lines xt^{\lambda}=const that focus around the origin. The existence of such solutions explores the invariance of the system (1) under rescalings and we look for profiles with U(\xi), \Gamma(\xi), \Theta(\xi) even functions and V(\xi) odd function.
We further demand that these profiles are localizing. We will call a selfsimilar function
f(t,x)=t^{b}F(xt^{\lambda})\,,\quad\mbox{with $F(\xi)=F(\xi)$ and $\lambda>0$}  (5) 
localizing if it has the asymptotic behavior
F(\xi)={{O}}(\xi^{p})\quad\mbox{ as $\xi\to\infty$ }  (6) 
and satisfies that p<0 when b>0 while p>0 when b<0. Under this definition, when f(t,0) grows then f(t,x) grows at a slower rate when x\neq 0, while when f(t,0) decays then f(t,x) decays at a slower rate at x\neq 0. We will call a selfsimilar function with an oddprofile F(\xi)=F(\xi) localizing when its derivative f_{x}(t,x) has the aforementioned behavior.
Applying the ansatz (4) leads to the system of singular ordinary differential equations
\displaystyle V^{\prime}(\xi)  \displaystyle=U(\xi),  (7)  
\displaystyle\Sigma^{\prime}(\xi)  \displaystyle=bV(\xi)+\lambda\xi U(\xi),  
\displaystyle a\Gamma(\xi)+\lambda\xi\Gamma^{\prime}(\xi)  \displaystyle=U(\xi),  
\displaystyle c\Theta(\xi)+\lambda\xi\Theta^{\prime}(\xi)  \displaystyle=\Sigma(\xi)U(\xi),  
\displaystyle\Sigma(\xi)  \displaystyle=\Theta(\xi)^{\alpha}\Gamma(\xi)^{m}U(\xi)^{n},  
\displaystyle\Gamma(0)=\Gamma_{0}>0,\quad U(0)  \displaystyle=U_{0}>0,\quad\text{$\xi\in[0,\infty)$}. 
This is viewed as a system of singular ordinary differential equations for (V,\Sigma,\Gamma,\Theta) with U defined by inverting (7){}_{5}. With the objective to compare these selfsimilar profiles to the fundamental solution of the heat or porous media equation, we look for solutions so that (\Gamma,\Theta) and U is even; in turn implying that \Sigma is even while V is odd (see section 3 for details). We impose
V(0)=U^{\prime}(0)=\Gamma^{\prime}(0)=\Sigma^{\prime}(0)=\Theta^{\prime}(0)=0  (8) 
so that the symmetric extensions are smooth selfsimilar profiles. The main result of this article is the construction of profiles solving (7), (8). It turns out that the induced selfsimilar solutions (4) exhibit localizing in space behavior as time evolves, see sections 8 and 9.
The idea of constructing selfsimilar localizing solutions for problems of shear band formation is introduced in [16] for the system
\displaystyle v_{t}  \displaystyle=(e^{\alpha\theta}v_{x}^{n})_{x}  (9)  
\displaystyle\theta_{t}  \displaystyle=e^{\alpha\theta}v_{x}^{n+1} 
modeling a nonNewtonian fluid with temperature dependent viscosity. Due to special properties of (9), the construction of selfsimilar solutions is reduced to finding a heteroclinic orbit for a planar system of autonomous differential equations, which is achieved through phase space analysis. A second step is taken in [22] where (1) is studied for parameters \alpha=0 and m<0 when the system simplifies to a system of two conservation laws. The problem is reduced using geometric singular perturbation theory to the flow of a planar dynamical system in a twodimensional invariant manifold. The present work addresses the full model (1) and due to the higher dimensionality of the problem a more elaborate version of geometric singular perturbation theory is needed.
The selfsimilar localizing solutions emerge as the combined outcome of Hadamard instability (that characterizes the system (1) for n=0 in the regime \alpha+m<0) and the regularizing effect of momentum diffusion when n>0. This feature can be clearly seen in the linearized analysis of uniform shearing solutions for the simplified model (9) which indicates that the combined effect of the two mechanisms amounts to Turing instability, see [16]. Moreover, existing linearized and nonlinear stability analyses that are available for special instances of (1) and are outlined in section 2 corroborate this point.
The article is organized as follows: Sections 3 and 4 deal with the formulation of the problem leading to (7), (8). The system (7) is singular (at \xi=0) and nonautonomous and it does not fit under a general existence theory. The singularity can be resolved and (7) is desingularized using again the scaleinvariance properties. Furthermore, upon introducing a series of nonlinear transformations, the construction of profiles for (7), (8) is accomplished by the construction of a heteroclinic orbit for the fourdimensional dynamical system for (p,q,r,s)
\displaystyle\dot{p}  \displaystyle=p\Big{(}\frac{1}{\lambda}(ra)+2\lambda prq\Big{)},  (S)  
\displaystyle\dot{q}  \displaystyle=q\Big{(}1\lambda prq\Big{)}+bpr,  
\displaystyle n\dot{r}  \displaystyle=r\Big{(}\frac{\alphamn}{\lambda(1+\alpha)}(ra)+\lambda pr+q+% \frac{\alpha}{\lambda}r\big{(}s\frac{1+m+n}{1+\alpha}\big{)}+\frac{n\alpha}{% \lambda(1+\alpha)}\Big{)},  
\displaystyle\dot{s}  \displaystyle=s\Big{(}\frac{\alphamn}{\lambda(1+\alpha)}(ra)+\lambda pr+q% \frac{1}{\lambda}r\big{(}s\frac{1+m+n}{1+\alpha}\big{)}\frac{n}{\lambda(1+% \alpha)}\Big{)}, 
parametrized by (\lambda,\alpha,m,n). The initial conditions are transmitted to asymptotic conditions for the heteroclinic as \eta(=\log\xi)\to\infty while the behavior as \eta\to\infty will capture the asymptotic behavior of the profiles.
The existence of solutions to (7), (8) is achieved in sections 3  7. Their construction is reduced to obtaining a heteroclinic orbit for (S) with prescribed asymptotic behavior as \eta\to\infty. At the end of section 3, the reader will find an outline on how the construction of the profiles is reduced to obtaining a heteroclinic orbit for (S). The existence of a heteroclinic orbit for (S) (with prescribed asymptotic behavior) is obtained in Theorem 1 using the geometric theory of singular perturbations [9, 10, 11, 12, 20, 21], exploiting the smallness of the parameter n. Section 7 contains the main part of the proof, motivated by the geometric arguments of [25] and adapted to the present system through somewhat cumbersome computations detailed in sections 7.2 and 7.3. The proof is based on a more elaborate argument than the simple invariant manifold argument for obtaining the corresponding result for the simplified model in [22]. In the present case, the finer structure inside the manifold is needed along with the persistence of the unstable and stable manifolds, see section 7.3.
The constructed selfsimilar solutions depend on two parameters (U_{0},\Gamma_{0}) describing the initial nonuniformity; the rate of localization \lambda is determined from (U_{0},\Gamma_{0}) via (48). Due to the construction necessities the rate has to obey the bound (48). The solutions (4) provide an example of instabilty resulting in localization. Their localizing behavior is investigated in section 8, see Proposition 8.1 and section 8.2. In section 9, the heteroclinic orbit is computed numerically using the standardized continuation software AUTO, [6, 7, 8], what leads to graphs of the profiles and the corresponding localizing solutions for various examples of material parameters.
To our knowledge, the localizing selfsimilar solutions are the first instance of depicting localizing behavior for a sufficiently broad model (1) that embodies the basic shear band formation mechanism proposed by Zener and Hollomon [31] and Clifton [2], and encompasses all the contributing factors of thermal softening, strain hardening and strainrate hardening. They complement [27], where shear bands are induced by energy supplied via the boundary. Some of the key predictions of stresscollapse are common, but the present result has the conceptual advantage to capture the emergence of localization as the combined result of Hadamard instability with small viscosity effects. It would be very interesting to study the stability of the solutions that are constructed here; this appears a challenging problem.
A preliminary report of these results, concerning the case with no strain hardening (m=0), has been presented in the Proceedings article [19].
2 Description of the shear band formation problem
The formation of shear bands [4, 31] is a phenomenon occuring during high strainrate plastic deformations of certain steels and other metal alloys. Instead of distributing evenly across the loaded region, the shear strain concentrates in a narrow band with a concurrent elevation of the temperature in the interior of the band, [31, 4, 14]. Shear bands are often precursors to rupture and their study has attracted considerable attention including experimental works [4, 14], mechanical modeling and linearized analysis studies (e.g. [3, 13, 23, 30] and references therein) and nonlinear analysis investigations [5, 27, 1].
2.1 Modeling shear bands
Shear bands appear and propagate as one dimensional structures (up to interaction times), and many investigations focus on the study of onedimensional, simple shear. A specimen located in the xyplain undergoes shear motion in the ydirection. The motion is described by the (plastic) shear strain \gamma(t,x), the strain rate u(t,x)=\gamma_{t}(t,x), the velocity v(t,x) in the shear direction, the temperature \theta(t,x) and the shear stress \sigma(t,x) all defined in (t,x)\in\mathbb{R}^{+}\times\mathbb{R}. It is described by the equations
\displaystyle\gamma_{t}  \displaystyle=v_{x}  (10)  
\displaystyle v_{t}  \displaystyle=\sigma_{x}  
\displaystyle\theta_{t}  \displaystyle=\kappa\theta_{xx}+\sigma v_{x}, 
which stand respectively for the kinematic compatibility equation, the balance of momentum and the balance of energy equation. Here, the elastic effects are neglected and all strain is considered to be plastic, and a Fourier heat conduction is considered with \kappa the thermal diffusivity.
Under shearing most materials deform in a uniform fashion until they break. By contrast, in high strainrate deformations of certain steels, it is observed that nonuniformities develop in the plastic strain and localize in a narrow region, called shear band; see Fig. 1 for a caricature of shear band forming. Shear bands correspond to material instabilities and are usually observed in the interior of specimens. Typically, the maximum temperature is measured in the interior of the band [4].
It was recognized by Zener and Hollomon [31] that the high deformation speed has two effects: First, an increase in the deformation speed changes the deformation conditions from isothermal to nearly adiabatic. Under such conditions the combined effect of thermal softening and strain hardening tends to produce net softening response. (Indeed, experimental observations of shear bands are typically associated with strain softening response – past a critical strain – of the measured stressstrain curve [3].) Second, strain rate has an effect per se, and needs to be included in the constitutive modeling.
Both effects are captured by modeling shear band formation via constitutive models within the framework of thermoviscoplasticity:
\sigma=f(\theta,\gamma,\gamma_{t})\quad\mbox{where}\quad f_{p}(\theta,\gamma,p% )>0\,.  (11) 
The constitutive relation (11) may be viewed as a yield surface or, upon inverting it, as a plastic flow rule. This suggests the terminology: the material exhibits thermal softening at state variables (\theta,\gamma,p) where f_{\theta}(\theta,\gamma,p)<0, strain hardening at state variables where f_{\gamma}(\theta,\gamma,p)>0, and strain softening when f_{\gamma}(\theta,\gamma,p)<0. The slopes of f – f_{\theta}, f_{\gamma} or f_{p} – measure respectively the degree of thermal softening, strain hardening (or softening) and strainrate sensitivity, respectively. The difficulty of performing high strainrate experiments causes uncertainty as to the specific form of the constitutive form of the stress. Here, we will use two constitutive laws to describe the stress \sigma:
\displaystyle\sigma=\theta^{\alpha}\gamma^{m}\gamma_{t}^{n},  \displaystyle\text{ power law },  (12)  
\displaystyle\sigma=e^{\alpha\theta}v_{x}^{n},  \displaystyle\text{ exponential law}\,.  (13) 
The power law (12) characterizes the response of the material. The parameter \alpha>0 measures the degree of thermal softening, m>0 measures the degree of strain hardening (or m<0 in case of a softening plastic flow), while n>0 measures strainrate hardening and is typically small, n\ll 1, [3, 2]. It is an empirical law and the parameters are determined by fitting experimental data.
We summarize the equations describing the model. For the power law the resulting system reads
\displaystyle v_{t}=\sigma_{x},  (14)  
\displaystyle\theta_{t}=\kappa\theta_{xx}+\sigma\gamma_{t},  
\displaystyle\gamma_{t}=v_{x},  
\displaystyle\sigma=\theta^{\alpha}\gamma^{m}\gamma_{t}^{n}\,. 
The system (14) captures the simplest mechanism proposed for shear localization in highspeed deformations of metals [31, 2], and an (isothermal) variant appears in early studies of necking [15]. Very often attention is restricted to the adiabatic model \kappa=0 which is appropriate for the initial development of shear bands under very fast deformations.
The exponential law does not exhibit any strain hardening and thus (10) decouples and leads to the simplified system
\displaystyle v_{t}=\sigma_{x},  (15)  
\displaystyle\theta_{t}=\kappa\theta_{xx}+\sigma v_{x}  
\displaystyle\sigma=e^{\alpha\theta}v_{x}^{n}\,. 
The exponential law can be interpreted as a temperature dependent nonNewtonian fluid and is exactly (9) for adiabatic deformations (\kappa=0).
2.2 Uniform shearing solutions
In the study of shear bands, a special problem is often considered where an infinite slab of material is sheared by prescribed constant velocity V=1 at the upper plate while the lower plate is held fixed. This is described by setting the plates at x=0,1 and imposing prescribed (normalized) velocities v(t,0)=0, v(t,1)=1, respectively. The plates are thermally insulated: \theta_{x}(t,0)=0, \theta_{x}(t,1)=0. For the heat flux Q one either uses the adiabatic assumption Q=0 (equivalently \kappa=0) or alternatively a Fourier law, Q=\kappa\theta_{x} with thermal diffusivity parameter \kappa. Imposing adiabatic conditions projects the belief that, at high strain rates, heat diffusion operates at a slower time scale than the timescale of the development of a shear band. It appears a plausible assumption for the shear band initiation process, but not necessarily for the evolution of a developed band, due to the high temperature differences involved.
The model (14) admits a special class of solutions describing uniform shearing: They emanate from spatially uniform initial data \gamma_{0} and \theta_{0}, and are obtained by the ansatz \gamma_{s}(t)=t+\gamma_{0} and v_{s}(x)=x for the strain and velocity respectively. They are obtained upon solving the ordinary differential equation
\frac{d\theta_{s}}{dt}=\sigma_{s}=\theta_{s}^{\alpha}(t+\gamma_{0})^{m}\,,% \quad\theta_{s}(0)=\theta_{0}\,,  (16) 
and read
\displaystyle v_{s}(x)  \displaystyle=x\,,\quad\gamma_{s}(t)=t+\gamma_{0},  (17)  
\displaystyle\theta_{s}(t)  \displaystyle=\left(\tfrac{1+\alpha}{1+m}\right)^{\frac{1}{1+\alpha}}(t+\gamma% _{0})^{\frac{1+m}{1+\alpha}}\left(1+\tfrac{1+m}{1+\alpha}\big{(}\theta_{0}^{1+% \alpha}\tfrac{1+\alpha}{1+m}\gamma_{0}^{1+m}\big{)}\tfrac{1}{(t+\gamma_{0})^{% m+1}}\right)^{\frac{1}{1+\alpha}}  
\displaystyle\sigma_{s}(t)  \displaystyle=\left(\tfrac{1+\alpha}{1+m}\right)^{\frac{\alpha}{1+\alpha}}(t+% \gamma_{0})^{\frac{\alpha+m}{1+\alpha}}\left(1+\tfrac{1+m}{1+\alpha}\big{(}% \theta_{0}^{1+\alpha}\tfrac{1+\alpha}{1+m}\gamma_{0}^{1+m}\big{)}\tfrac{1}{(t% +\gamma_{0})^{m+1}}\right)^{\frac{\alpha}{1+\alpha}}. 
Equation (17){}_{3} describes the stressstrain curve \sigma_{s} versus \gamma_{s} for uniform shear. The stressstrain curve is increasing when \alpha<m but it is decreasing for large times when \alpha>m. Here, we are interested in the regime \alpha>m where thermal softening dominates strain hardening and produces net softening.
2.3 On the stability of the uniform shearing solution
The system (1) for n=0 is a firstorder system. When \alpha>m, the initial value problem has two purely imaginary eigenvalues in a regime of strain beyond the maximum of the stressstrain curve (see Appendix A). Accordingly, the linearized system (for n=0) around the uniform shearing solution (17) exhibits Hadamard instability, see Appendix A.
The stability of the uniform shear solution for n>0 has been the objective of many investigations. Since (17) is time dependent this leads to investigations of nonautonomous systems. A natural way to define stability is to consider
\displaystyle\gamma_{s}^{*}(t)=t+\gamma_{0},\quad\theta_{s}^{*}(t)  \displaystyle=\left(\tfrac{1+\alpha}{1+m}\right)^{\frac{1}{1+\alpha}}(t+\gamma% _{0})^{\frac{1+m}{1+\alpha}}\,,  (18) 
the functions capturing the growth of the uniform shearing solution, and to study the relative perturbations
u(t,x)=v_{x}(t,x)\,,\quad\hat{\Gamma}(t,x)=\frac{\gamma(t,x)}{\gamma_{s}^{*}(t% )}\,,\quad\hat{\Theta}(t,x)=\frac{\theta(t,x)}{\theta_{s}^{*}(t)}.  (19) 

The uniform shear solution is asymptotically stable when the solution emanating from small perturbations of (17) satisfies that (u,\hat{\Gamma},\hat{\Theta})\to(1,1,1) as time goes to infinity.

The uniform shear solution is unstable if for small perturbations of (17) the relative perturbations (u,\hat{\Gamma},\hat{\Theta}) drift away from (1,1,1) as time increases.
This notion of stability is used in nonlinear stability studies of shear bands [5, 26] as well as in linearized stability analyses by Molinari and Clifton [23, 13] who coined the name stability analysis of relative perturbations. The problem of stability is presently resolved only for the special cases m=0 or \alpha=0 for (1); these are cases that the system decouples and reduces to simpler models:
Understanding of the nature of the instability is offered in [16] for the model (9), which has the special property that both the nonlinear and the linearized analysis of relative perturbations is reduced to studying autonomous systems. In particular, linearized stability (or instability) can be accessed via analyzing Fourier modes; see [16]. For n=0, the linearized stability analysis predicts exponential growth of the high frequency modes, leading to what is usually termed as Hadamard instability. By contrast, when n>0 the linear modes are still unstable and their growth rates are increasing with frequency but they are uniformly bounded by a bound independent of the frequency. The behavior of the linearized system around the uniform shearing solution for the full system (1) is at present open; the conjecture is that it has the same structure as described above for relative perturbations of (9) when n>0 is small, and it is stable past a certain threshold. This is corroborated by linearized analysis for the special case (1) with \alpha=0, m=1, n<1, which again leads to the study of autonomous systems for relative perturbations, [18].
2.4 The nonlinear regime
In the unstable parameter regime, at the initial stage unstable modes start to grow and this process can be captured by the linearized problem. The second stage of localization lies within the realm of nonlinear analysis. The question arises how the high frequency oscillations resulting from Hadamard instability interact with the nonlinearity and the viscosity to form a coherent structure. An asymptotic criterion accounting for the nonlinear aspects of localization is derived in [17]. Based on ideas from the theory of relaxation system and the Chapman Enskog expansion, an effective equation is derived for the nonlinear dynamics (1). It predicts stability in the regime \alpha+m+n>0 and instability in the regime \alpha+m+n<0, see [17].
Insight on how coherent structures form can be offered by investigating selfsimilar solutions (4). It is customary in studies of parabolic systems (like (3)) to investigate diffusing selfsimilar solutions corresponding to the parameter selection \lambda<0. By contrast, selfsimilar solutions with \lambda>0 tend to propagate information along the lines t^{\lambda}x=const. and thus to localize around the point x=0. Selfsimilar localizing solutions were established in [16] for the model (9) using a phaseplane analysis for the resulting twodimensional system. They will be pursued also here for the power law (1).
3 Selfsimilar solutions
We consider the system (1) (or the system (3)) in the domain x\in\mathbb{R}, t>0 and note that both systems are invariant under a family of scaling transformations: if (\gamma,u,v,\theta,\sigma) satisfy (1), with u, \sigma connected via (2), then for any \lambda\in\mathbb{R} and \rho>0 the rescaled functions (\gamma_{\rho},u_{\rho},v_{\rho},\theta_{\rho},\sigma_{\rho}) defined by
\displaystyle\gamma_{\rho}(t,x)  \displaystyle=\rho^{a}\gamma(\rho^{1}t,\rho^{\lambda}x),  \displaystyle v_{\rho}(t,x)  \displaystyle=\rho^{b}v(\rho^{1}t,\rho^{\lambda}x),  (20)  
\displaystyle\theta_{\rho}(t,x)  \displaystyle=\rho^{c}\theta(\rho^{1}t,\rho^{\lambda}x),  \displaystyle\sigma_{\rho}(t,x)  \displaystyle=\rho^{d}\sigma(\rho^{1}t,\rho^{\lambda}x),  
\displaystyle u_{\rho}(t,x)  \displaystyle=\rho^{b+\lambda}\gamma(\rho^{1}t,\rho^{\lambda}x) 
also satisfies (1), provided
\displaystyle a  \displaystyle:=a_{0}+a_{1}\lambda=\frac{2+2\alphan}{D}+\frac{2(1+\alpha)}{D}\lambda,  \displaystyle b  \displaystyle:=b_{0}+b_{1}\lambda=\frac{1+m}{D}+\frac{1+m+n}{D}\lambda,  (21)  
\displaystyle c  \displaystyle:=c_{0}+c_{1}\lambda=\frac{2(1+m)}{D}+\frac{2(1+m+n)}{D}\lambda,  \displaystyle d  \displaystyle:=d_{0}+d_{1}\lambda=\frac{2\alpha+2m+n}{D}+\frac{2(\alpha+m+n)% }{D}\lambda, 
and
D=1+2\alphamn\,.  (22) 
The same scaling trasformation leaves invariant solutions of (3). We note there are two independent scaling parameters in (20), \rho and \lambda, while the remaining parameters are determined by the relations (21), (22). Throughout this work, the material parameters (\alpha,m,n) will be restricted to the range
\displaystyle\alpha>0  \displaystyle\text{(thermal softening)},  (23)  
\displaystyle m>1  \displaystyle\text{(strain softening/hardening)},  
\displaystyle n>0  \displaystyle\text{(strain rate sensitivity)},  
\displaystyle\alpha+m+n<0  \displaystyle\text{(unstable regime)}. 
Observe that (23){}_{4} implies that \alpha+m<0 and thus we are in the regime of net softening, where the associated hyperbolic system with n=0 loses hyperbolicity, see Appendix A. Moreover, D>1+\alpha>1 while 1+\alphan>1+m>0.
Solutions of (1) or (3) that are selfsimilar with respect to the scaling transformation (20), (21), (22) have the form
\displaystyle\gamma(t,x)  \displaystyle=t^{a}\Gamma(t^{\lambda}x),  \displaystyle v(t,x)  \displaystyle=t^{b}V(t^{\lambda}x),  \displaystyle\theta(t,x)  \displaystyle=t^{c}\Theta(t^{\lambda}x),  (24)  
\displaystyle\sigma(t,x)  \displaystyle=t^{d}\Sigma(t^{\lambda}x),  \displaystyle u(t,x)  \displaystyle=t^{b+\lambda}U(t^{\lambda}x)\,, 
and depend on one parameter, \lambda. In the sequel, we are interested in constructing solutions (24) defined in the domain x\in\mathbb{R}, t>0 for values of the parameter \lambda>0.
To motivate the role of selfsimilar solutions with \lambda>0 and some forthcoming selections, recall that the fundamental solution of the heat equation u_{t}=u_{xx} is of selfsimilar form
u(t,x)=\frac{1}{\sqrt{t}}U(\frac{x}{\sqrt{t}})\,. 
Moreover, power nonlinear parabolic diffusion equations (such as the porous media) admit selfsimilar solutions which correspond to values \lambda<0 and capture the effect of diffusion. We are interested here to investigate whether the couplings with the remaining equations in (3) can lead to the opposite behavior, of localization, and we seek existence of selfsimilar solutions with the parameter in the range \lambda>0. Note that profiles of the form (24) with \lambda>0 are constant on lines \xi=t^{\lambda}x and are thus expected to localize in space as time evolves. In order to compare the solutions we intend to construct for \lambda>0, with the existing selfsimilar solutions of nonlinear parabolic equations, we seek solutions where u(t,x) admits a maximum located at x=0 for all times. This imposes for selfsimilar solutions that
U^{\prime}(0)=0  (25) 
and as we will see this induces some symmetry properties to solutions.
Introducing the ansatz (24) into (1) gives a system of ordinary differential equations
\displaystyle V^{\prime}  \displaystyle=U  (26)  
\displaystyle\Sigma^{\prime}  \displaystyle=bV+\lambda\xi U  
\displaystyle c\Theta+\lambda\xi\Theta^{\prime}  \displaystyle=\Sigma U  
\displaystyle a\Gamma+\lambda\xi\Gamma^{\prime}  \displaystyle=U 
together with an algebraic equation,
\Sigma=\Theta^{\alpha}\Gamma^{m}U^{n}\,,  (27) 
obtained from (2). This is viewed as a first order system for the variable (V,\Sigma,\Theta,\Gamma)(\xi) with U(\xi) determined by inverting (27). In principle, solutions of (26) will depend on five data inputs: the initial data (V_{0},\Sigma_{0},\Theta_{0},\Gamma_{0}) and the parameter \lambda. The system (26) is nonautonomous and singular at \xi=0, what imposes compatibility conditions.
For smooth initial data, solutions of (3) either blowup in finite time or they are as smooth as the initial data [27, Theorem 1] (in fact analytic for analytic initial data). Thus the selfsimilar solutions will be sought to be smooth with their maxima fixed at the origin. The latter can be always achieved due to the translation invariance of (3). Next, we discuss the conditions imposed by these requirements: The initial conditions
V(0)=V_{0}\in\mathbb{R}\,,\quad\Sigma(0)=\Sigma_{0}>0,\quad\Theta(0)=\Theta_{0% }>0\,,\quad\Gamma(0)=\Gamma_{0}>0\,,\quad U(0)=U_{0}>0  (28) 
are supplemented with (25). Since the solution is smooth the singularity at \xi=0 imposes two compatibility conditions on the data which together with (27) imply
a\Gamma_{0}=U_{0}\,,\quad c\Theta_{0}=\Sigma_{0}U_{0}\,,\quad\Sigma_{0}=\Theta% _{0}^{\alpha}\Gamma_{0}^{m}U_{0}^{n}  (29) 
with a,b,c given by (21). The condition U^{\prime}(0)=0 together with the smoothness of the solution yields upon differentiating (26) and (27)
\displaystyle(a+\lambda)\Gamma^{\prime}(0)  \displaystyle=U^{\prime}(0)=0  (30)  
\displaystyle(c+\lambda)\Theta^{\prime}(0)  \displaystyle=\Sigma^{\prime}(0)U(0)+\Sigma(0)U^{\prime}(0)  
\displaystyle=\alpha\frac{\Sigma_{0}U_{0}}{\Theta_{0}}\Theta^{\prime}(0)+m% \frac{\Sigma_{0}U_{0}}{\Gamma_{0}}\Gamma^{\prime}(0). 
By (21) and (23), for \lambda>0 we have a>0, c>0, hence
\Gamma^{\prime}(0)=0\,,\quad\Theta^{\prime}(0)=0\,.  (31) 
Again by (25), (27), (26){}_{2} and b>0 ,
\Sigma^{\prime}(0)=0\,,\quad V(0)=V_{0}=0\,.  (32) 
Finally, since (26) is invariant under the change of variables
\xi\to\xi\,,\quad V\toV\,,\quad U\to U\,,\quad\Theta\to\Theta\,,\quad\Gamma\to\Gamma 
it admits solutions such that U, \Theta, \Gamma and \Sigma are even functions of \xi, while V is an odd function of \xi.
In summary, we proceed as follows: We first construct a solution (V(\xi),\Sigma(\xi),\Theta(\xi),\Gamma(\xi)) of (26) defined for \xi\in[0,\infty) and set U(\xi) by (27). The solution will be sought subject to the data
U^{\prime}(0)=\Gamma^{\prime}(0)=\Sigma^{\prime}(0)=\Theta^{\prime}(0)=0  (33) 
V(0)=0,\quad U(0)=U_{0}>0\,,\quad\Gamma(0)=\Gamma_{0}>0\,,\quad\Theta(0)=% \Theta_{0}>0\,,\quad\Sigma(0)=\Sigma_{0}>0\,.  (34) 
satisfying the compatibility conditions (29) for some \lambda>0. It is not apriori clear that this problem is not overdetermined and we give a detailed analysis of this point in Section 6. The constructed solution is then extended on (\infty,0] by setting
V(\xi)=V(\xi)\,,\quad U(\xi)=U(\xi)\,,\quad\Theta(\xi)=\Theta(\xi)\,,\quad% \Gamma(\xi)=\Gamma(\xi)\,, 
that is we use an odd extension for V and even extensions for U, \Theta, \Gamma and \Sigma. Given the material constants (\alpha,m,n) there are two independent parameters in the problem, which may be viewed as \Gamma_{0}, U_{0} and \lambda>0 subject to the constraint
\frac{U_{0}}{\Gamma_{0}}=\frac{2+2\alphan}{D}+\frac{2(1+\alpha)}{D}\lambda\,.  (35) 
The remaining constants \Theta_{0} and \Sigma_{0} are determined via (29).
The profiles are constructed in the forthcoming sections 47. Then in section 8 we check that the profiles are localizing according to the definition (5)(6) in the Introduction. This is based on the asymptotic behavior of the constructed profiles \xi\to\infty, established in Proposition 8.1.
Note that the uniform shearing solution is achieved as a selfsimilar profile for \lambda=\frac{1+m}{2(1+\alpha)}<0 and
\Gamma(\xi)=U(\xi)=U_{0},\quad V(\xi)=U_{0}\xi,\quad\Theta(\xi)=\Big{(}\frac{1% +\alpha}{1+m}U_{0}^{1+m+n}\Big{)}^{\frac{1}{1+\alpha}},\quad\Sigma(\xi)=\Big{(% }\frac{1+\alpha}{1+m}\Big{)}^{\frac{\alpha}{1+\alpha}}U_{0}^{\frac{\alpha+m+% n}{1+\alpha}}. 
The uniform shear should be contrasted to the solutions that are constructed here which exhibit localizing behavior: the growth of the strain is superlinear at the origin and the profiles of the solution (at fixed times) localize as time proceeds, see Section 8.
We give a short roadmap of how we proceed to construct the solution of (26), (33), (34) and determine its properties.
As an outcome of this construction, it turns out there is a two parameter family of solutions depending on the data U_{0} and \Gamma_{0} with the rate \lambda determined via (35). The dynamic stability of the solutions is a challenging open problem.
4 Reduction to the construction of a heteroclinic orbit
The goal of this section is to derive an equivalent system (S) to (26) that is autonomous and to turn the problem of constructing profiles for (26) to the construction of a heteroclinic orbit for (S). We employ techniques from [16] and [22]. The novelty of the present analysis lies in the higher dimensionality of the resulting system especially with regard to the construction of the heteroclinic orbit.
4.1 Desingularization
We regard (26) as a boundaryvalue problem in the right halfline \xi\in[0,\infty) subject to the boundary conditions (33) and proceed to desingularize it. The system (26) is itself scale invariant: Given a solution \big{(}\Gamma(\xi),V(\xi),\Theta(\xi),\Sigma(\xi),U(\xi)\big{)} the rescaled function \big{(}\Gamma_{\rho}(\xi),V_{\rho}(\xi),\Theta_{\rho}(\xi),\Sigma_{\rho}(\xi),% U_{\rho}(\xi)\big{)} defined by
\displaystyle\Gamma_{\rho}(\xi)  \displaystyle=\rho^{a_{1}}\Gamma(\rho\xi),  \displaystyle V_{\rho}(\xi)  \displaystyle=\rho^{b_{1}}V(\rho\xi),\quad\Theta_{\rho}(\xi)=\rho^{c_{1}}% \Theta(\rho\xi),  (36)  
\displaystyle\Sigma_{\rho}(\xi)  \displaystyle=\rho^{d_{1}}\Sigma(\rho\xi),  \displaystyle U_{\rho}(\xi)  \displaystyle=\rho^{b_{1}+1}U(\rho\xi)=\rho^{a_{1}}U(\rho\xi) 
is again a solution. The class of functions that remain invariant under this scaling transformation is
\big{(}\Gamma(\xi),V(\xi),\Theta(\xi),\Sigma(\xi),U(\xi)\big{)}=\big{(}A\xi^{% a_{1}},B\xi^{b_{1}},C\xi^{c_{1}},D\xi^{d_{1}},E\xi^{a_{1}}\big{)} 
where A,B,C,D,E constants. Such a function is singular at \xi=0 and fails to satisfy (33). Nevertheless, it suggests the change of variables
\displaystyle{\bar{\gamma}}(\xi)  \displaystyle=\xi^{a_{1}}\Gamma(\xi),  \displaystyle{\bar{v}}(\xi)  \displaystyle=\xi^{b_{1}}V(\xi),  \displaystyle{\bar{\theta}}(\xi)  \displaystyle=\xi^{c_{1}}\Theta(\xi),  \displaystyle{\bar{\sigma}}(\xi)  \displaystyle=\xi^{d_{1}}\Sigma(\xi),  \displaystyle{\bar{u}}(\xi)  \displaystyle=\xi^{b_{1}+1}U(\xi),  (37) 
with a_{1},b_{1},c_{1} and d_{1} as in (21), in order to desingularize the problem. After some cumbersome, but straightforward calculation, we find that ({\bar{\gamma}},{\bar{v}},{\bar{\theta}},{\bar{\sigma}}) satisfies
\displaystyle a_{0}{\bar{\gamma}}+\lambda\xi{\bar{\gamma}}^{\prime}  \displaystyle={\bar{u}},  (38)  
\displaystyle b_{0}{\bar{v}}+\lambda\xi{\bar{v}}^{\prime}  \displaystyle=d_{1}{\bar{\sigma}}+\xi{\bar{\sigma}}^{\prime},  
\displaystyle c_{0}{\bar{\theta}}+\lambda\xi{\bar{\theta}}^{\prime}  \displaystyle={\bar{\sigma}}{\bar{u}},  
\displaystyleb_{1}{\bar{v}}+\xi{\bar{v}}^{\prime}  \displaystyle={\bar{u}}, 
and {\bar{u}} is defined by
{\tilde{\sigma}}={\bar{\theta}}^{\alpha}{\bar{\gamma}}^{m}{\bar{u}}^{n}. 
Next, introduce a new independent variable \eta=\log\xi and define ({\tilde{\gamma}},{\tilde{v}},{\tilde{\theta}},{\tilde{\sigma}},{\tilde{u}}) by
\displaystyle{\tilde{\gamma}}(\log\xi)  \displaystyle={\bar{\gamma}}(\xi),  \displaystyle{\tilde{v}}(\log\xi)  \displaystyle={\bar{v}}(\xi),  \displaystyle{\tilde{\theta}}(\log\xi)  \displaystyle={\bar{\theta}}(\xi),  (39)  
\displaystyle{\tilde{\sigma}}(\log\xi)  \displaystyle={\bar{\sigma}}(\xi),  \displaystyle{\tilde{u}}(\log\xi)  \displaystyle={\bar{u}}(\xi). 
Noticing that \frac{d}{d\eta}{\tilde{\gamma}}(\eta)=\xi\frac{d}{d\xi}{\bar{\gamma}}(\xi), we obtain an autonomous system
\displaystyle a_{0}{\tilde{\gamma}}+\lambda{\dot{\tilde{\gamma}}}  \displaystyle={\tilde{u}},  (40)  
\displaystyle b_{0}{\tilde{v}}+\lambda{\dot{\tilde{v}}}  \displaystyle=d_{1}{\tilde{\sigma}}+{\dot{\tilde{\sigma}}},  
\displaystyle c_{0}{\tilde{\theta}}+\lambda{\dot{\tilde{\theta}}}  \displaystyle={\tilde{\sigma}}{\tilde{u}},  
\displaystyleb_{1}{\tilde{v}}+{\dot{\tilde{v}}}  \displaystyle={\tilde{u}}, 
where the notation \dot{f}=\frac{df}{d\eta} is used, and {\tilde{u}} is defined by
{\tilde{\sigma}}={\tilde{\theta}}^{\alpha}{\tilde{\gamma}}^{m}{\tilde{u}}^{n}\,. 
The system (40) is autonomous and one might attempt to consider its equilibria. However, it is easy to conclude that we cannot expect a heteroclinic that tends to equilibria of (40). Indeed, suppose {\tilde{u}}\rightarrow{\tilde{u}}_{\infty}\geq 0 as \eta\rightarrow\infty. Then from the last equation in (40), we conclude that {\tilde{v}}\rightarrow\infty. This suggests to enlarge the scope and consider solutions that grow as polynomials (or faster) at infinities.
4.2 The (p,q,r,s)system derivation
Next, we attempt to come up with a new choice of variables that tend to equilibria as \eta\rightarrow\pm\infty and accommodate orbits that have power behavior at infinities. We rewrite (40) in the form
\displaystyle\frac{d}{d\eta}{(\ln{{\tilde{\gamma}}})}  \displaystyle=\tfrac{1}{\lambda}\big{(}a_{0}+\frac{{\tilde{u}}}{{\tilde{% \gamma}}}\big{)},  (41)  
\displaystyle\frac{d}{d\eta}{(\ln{{\tilde{v}}})}  \displaystyle=b_{1}+\frac{{\tilde{u}}}{{\tilde{v}}},  
\displaystyle\frac{d}{d\eta}{(\ln{{\tilde{\theta}}})}  \displaystyle=\tfrac{1}{\lambda}\big{(}c_{0}+\frac{{\tilde{\sigma}}{\tilde{u}% }}{{\tilde{\theta}}}\big{)},  
\displaystyle\frac{d}{d\eta}{(\ln{{\tilde{\sigma}}})}  \displaystyle=d_{1}+b\frac{{\tilde{v}}}{{\tilde{\sigma}}}+\lambda\frac{{\tilde% {u}}}{{\tilde{\sigma}}} 
and view it as describing the evolution of ({\tilde{\gamma}},{\tilde{v}},{\tilde{\theta}},{\tilde{\sigma}}) with {\tilde{u}} determined by {\tilde{u}}=\left(\frac{{\tilde{\sigma}}}{{\tilde{\theta}}^{\alpha}{\tilde{% \gamma}}^{m}}\right)^{\frac{1}{n}}.
This leads us to define
\displaystyle p:=\frac{{\tilde{\gamma}}}{{\tilde{\sigma}}},\quad q:=b\frac{{% \tilde{v}}}{{\tilde{\sigma}}},\quad r=\left(\frac{{\tilde{\sigma}}}{{\tilde{% \theta}}^{\alpha}{\tilde{\gamma}}^{m+n}}\right)^{\frac{1}{n}}=\frac{{\tilde{u% }}}{{\tilde{\gamma}}},\quad s:=\frac{{\tilde{\sigma}}{\tilde{\gamma}}}{{\tilde% {\theta}}}\,.  (42) 
The transformation (p,q,r,s)\leftrightarrow({\tilde{\gamma}},{\tilde{v}},{\tilde{\theta}},{\tilde% {\sigma}}) is a bijection in the positive orthant with the inverse determined by
{\tilde{\gamma}}=p^{\frac{1+\alpha}{D}}s^{\frac{\alpha}{D}}r^{\frac{n}{D}}\,,% \quad{\tilde{\theta}}=p^{\frac{1+m+n}{D}}s^{\frac{m+n1}{D}}r^{\frac{2n}{D}}\,, 
and then
{\tilde{\sigma}}=\frac{1}{{\tilde{\gamma}}}p\,,\quad{\tilde{v}}=\frac{1}{b}{% \tilde{\sigma}}\,q=\frac{1}{b}\frac{pq}{{\tilde{\gamma}}}\,. 
\displaystyle\frac{\dot{p}}{p}  \displaystyle=\frac{{\dot{\tilde{\gamma}}}}{{\tilde{\gamma}}}\frac{{\dot{% \tilde{\sigma}}}}{{\tilde{\sigma}}}  \displaystyle=\left[\frac{1}{\lambda}\Big{(}\frac{{\tilde{u}}}{{\tilde{\gamma}% }}a_{0}\Big{)}\right]  \displaystyle\left[d_{1}+b\frac{{\tilde{v}}}{{\tilde{\sigma}}}+\lambda\frac{{% \tilde{u}}}{{\tilde{\gamma}}}\frac{{\tilde{\gamma}}}{{\tilde{\sigma}}}\right]  
\displaystyle\frac{\dot{q}}{q}  \displaystyle=\frac{{\dot{\tilde{v}}}}{{\tilde{v}}}\frac{{\dot{\tilde{\sigma}% }}}{{\tilde{\sigma}}}  \displaystyle=\left[b_{1}+\frac{{\tilde{u}}}{{\tilde{v}}}\right]  \displaystyle\left[d_{1}+b\frac{{\tilde{v}}}{{\tilde{\sigma}}}+\lambda\frac{{% \tilde{u}}}{{\tilde{\gamma}}}\frac{{\tilde{\gamma}}}{{\tilde{\sigma}}}\right]  
\displaystyle n\frac{\dot{r}}{r}  \displaystyle=(m+n)\frac{{\dot{\tilde{\gamma}}}}{{\tilde{\gamma}}}+\frac{{% \dot{\tilde{\sigma}}}}{{\tilde{\sigma}}}+\alpha\frac{{\dot{\tilde{\theta}}}}{{% \tilde{\theta}}}  \displaystyle=\left[\frac{(m+n)}{\lambda}\Big{(}\frac{{\tilde{u}}}{{\tilde{% \gamma}}}a_{0}\Big{)}\right]  \displaystyle+\left[d_{1}+b\frac{{\tilde{v}}}{{\tilde{\sigma}}}+\lambda\frac{{% \tilde{u}}}{{\tilde{\gamma}}}\frac{{\tilde{\gamma}}}{{\tilde{\sigma}}}\right]+% \left[\frac{\alpha}{\lambda}\Big{(}\frac{{\tilde{\sigma}}{\tilde{u}}}{{\tilde{% \theta}}}c_{0}\Big{)}\right]  
\displaystyle\frac{\dot{s}}{s}  \displaystyle=\frac{{\dot{\tilde{\gamma}}}}{{\tilde{\gamma}}}+\frac{{\dot{% \tilde{\sigma}}}}{{\tilde{\sigma}}}\frac{{\dot{\tilde{\theta}}}}{{\tilde{% \theta}}}  \displaystyle=\left[\frac{1}{\lambda}\Big{(}\frac{{\tilde{u}}}{{\tilde{\gamma}% }}a_{0}\Big{)}\right]  \displaystyle+\left[d_{1}+b\frac{{\tilde{v}}}{{\tilde{\sigma}}}+\lambda\frac{{% \tilde{u}}}{{\tilde{\gamma}}}\frac{{\tilde{\gamma}}}{{\tilde{\sigma}}}\right]% \left[\frac{1}{\lambda}\Big{(}\frac{{\tilde{\sigma}}{\tilde{u}}}{{\tilde{% \theta}}}c_{0}\Big{)}\right]. 
We note that
\displaystyle\frac{{\tilde{\sigma}}{\tilde{u}}}{{\tilde{\theta}}}=rs,\quad% \frac{{\tilde{u}}}{{\tilde{v}}}=\frac{bpr}{q},\quad\frac{{\tilde{u}}}{{\tilde{% \sigma}}}=pr, 
and using (21) and (22), after a cumbersome but straightforward calculation, we derive the (p,q,r,s)system:
\displaystyle\dot{p}  \displaystyle=p\Big{(}\frac{1}{\lambda}(ra)+2\lambda prq\Big{)},  (S)  
\displaystyle\dot{q}  \displaystyle=q\Big{(}1\lambda prq\Big{)}+bpr,  
\displaystyle n\dot{r}  \displaystyle=r\Big{(}\frac{\alphamn}{\lambda(1+\alpha)}(ra)+\lambda pr+q+% \frac{\alpha}{\lambda}r\big{(}s\frac{1+m+n}{1+\alpha}\big{)}+\frac{n\alpha}{% \lambda(1+\alpha)}\Big{)},  
\displaystyle\dot{s}  \displaystyle=s\Big{(}\frac{\alphamn}{\lambda(1+\alpha)}(ra)+\lambda pr+q% \frac{1}{\lambda}r\big{(}s\frac{1+m+n}{1+\alpha}\big{)}\frac{n}{\lambda(1+% \alpha)}\Big{)}. 
In the sequel, we analyze (S) as an autonomous system: We begin with sorting its equilibria and analyzing their linear stability. Most importantly, (S) possesses the fastslow structure because of the small parameter n in the lefthandside of \eqref{eq:slow}_{3}; the dynamics of r can be distinctively faster than those of the other variables off the nullcline of r.
5 Equilibria and their linear stability
System (S) admits several equilibria listed in Appendix B. Our region of interest is the sector
\mathrsfs{P}=\{(p,q,r,s)\;\;p\geq 0,q\geq 0,r>0,s>0\}. 
That p,q\geq 0 comes from the requirement that {\tilde{\gamma}},{\tilde{v}},{\tilde{\sigma}}\geq 0. The reason we restrict to r>0, s>0 stems from mechanical considerations: If we transform back to the original variables, then we find
r(\eta)_{\eta=t^{\lambda}x}=t\partial_{t}\log\gamma(t,x),\quad r(\eta)s(\eta)% _{\eta=t^{\lambda}x}=t\partial_{t}\log\theta(t,x). 
Shear band initiation is related to conditions of loading where both the plastic strain and the temperature are increasing. This motivates to restrict to selfsimilar solutions taking values in the region r>0, s>0.
From the complete set of equilibria for (S) listed in Appendix B only two reside in the region r>0, s>0, namely
\displaystyle M_{0}  \displaystyle=(0,0,r_{0},s_{0}),  \displaystyle r_{0}  \displaystyle=a,  \displaystyle s_{0}  \displaystyle=\frac{1+m+n}{1+\alpha}\frac{n}{(1+\alpha)r_{0}},  
\displaystyle M_{1}  \displaystyle=(0,1,r_{1},s_{1}),  \displaystyle r_{1}  \displaystyle=a\frac{1+\alpha}{\alphamn}\lambda,  \displaystyle s_{1}  \displaystyle=\frac{1+m+n}{1+\alpha}\frac{n}{(1+\alpha)r_{1}}\,. 
\displaystyle a=a_{0}+a_{1}\lambda  \displaystyle=\frac{2+2\alphan}{D}+\frac{2+2\alpha}{D}\lambda  
\displaystyle D  \displaystyle=1+2\alphamn\,, 
and that the parameters (\alpha,m,n) take values in the range (23). As a consequence r_{0}>0, and a simple calculation shows that r_{0}s_{0}=\frac{2(1+m)}{D}+\frac{2(1+m+n)}{D}\lambda>0; hence, r_{0},s_{0}>0 and M_{0} resides in the region \mathrsfs{P}. By contrast, M_{1} can be out of the region r>0, s>0 if \lambda is large enough. Note that r_{1},s_{1}>0 only if \frac{1+m+n}{1+\alpha}r_{1}>\frac{n}{(1+\alpha)}. This reads
\frac{1+m+n}{1+\alpha}\Big{(}\frac{2+2\alphan}{D}\frac{(1+\alpha)(1+m+n)}{D(% \alphamn)}\lambda\Big{)}>\frac{n}{(1+\alpha)}, 
and thus M_{1} resides in the region \mathrsfs{P} only under the constraint
0<\lambda<\frac{2(\alphamn)}{1+m+n}\left(\frac{1+m}{1+m+n}\right).  (43) 
Henceforth, we restrict attention to rates \lambda satisfying (43).
We denote the four eigenvalues and four eigenvectors of the vector field linearized at M_{i}, i=0,1, by \mu_{ij} and X_{ij} with j=1,2,3,4.

M_{0} is a saddle; the matrix of the linearized vector field at M_{0} has three positive eigenvalues and one negative eigenvalue.
\mu_{01}=2,\quad\mu_{02}=1,\quad\mu_{03}=\mu_{0}^{+}={{O}}\Big{(}\frac{1}{n}% \Big{)}>0,\quad\mu_{04}=\mu_{0}^{}<0, (44) where \mu_{0}^{\pm} are respectively a positive and a negative solution of the quadratic equation
\Big{(}\mu\frac{r_{0}}{n}\Big{(}\frac{1s_{0}}{\lambda}\frac{n}{\lambda r_{0% }}\Big{)}\Big{)}\Big{(}\mu+\frac{s_{0}r_{0}}{\lambda}\Big{)}\frac{s_{0}r_{0}(% 1s_{0})(\alpha r_{0})}{n\lambda^{2}}=0. The leading orders of \mu_{0}^{\pm} are given by
\ \mu_{0}^{+}=\frac{1}{n}\tfrac{2(\alpham)(1+\alpha)(1+\lambda)}{\lambda(1+% \alpha)(1+2\alpham)}+{{O}}(1),\quad\mu_{0}^{}=\tfrac{2(1+m)(1+\alpha)(1+% \lambda)}{\lambda(1+2\alpham)}+{{O}}(n). Notice that one of the positive eigenvalue \mu_{03} is {{O}}(\frac{1}{n}), which indicates the separably fast dynamics along the direction X_{03}. We will make use of this structure later. The precise eigenvector components are presented in Appendix C, the directions of the eigenvectors are pointed out in Fig. 2 for n sufficiently small.

M_{1} is a saddle; the matrix of the linearized vector field at M_{1} has one positive eigenvalue and three negative eigenvalues.
\mu_{11}=\frac{1+m+n}{\alphamn},\quad\mu_{12}=1,\quad\mu_{13}=\mu_{1}^{+}=% {{O}}\Big{(}\frac{1}{n}\Big{)}>0,\quad\mu_{14}=\mu_{1}^{}<0, (45) where \mu_{1}^{\pm} is respectively a positive and a negative solution of the quadratic equation
\Big{(}\mu\frac{r_{1}}{n}\Big{(}\frac{1s_{1}}{\lambda}\frac{n}{\lambda r_{1% }}\Big{)}\Big{)}\Big{(}\mu+\frac{s_{1}r_{1}}{\lambda}\Big{)}\frac{s_{1}r_{1}(% 1s_{1})(\alpha r_{1})}{n\lambda^{2}}=0. The leading orders of \mu_{1}^{\pm} are given by
\displaystyle\mu_{1}^{+} \displaystyle=\frac{\alpham}{n\lambda(1+\alpha)}\Big{(}\frac{2(1+\alpha)(1+% \lambda)}{(1+2\alpham)}\frac{1+\alpha}{\alphamn}\lambda\Big{)}+{{O}}(1), \displaystyle\mu_{1}^{} \displaystyle=\frac{1+m}{\lambda}\Big{(}\frac{2(1+\alpha)(1+\lambda)}{(1+2% \alpham)}\frac{1+\alpha}{\alphamn}\lambda\Big{)}+{{O}}(n). Note that the positive eigenvalue \mu_{13} is {{O}}(\frac{1}{n}). In constrast to what happens at M_{0}, the eigenvalues of the linearized vector field at M_{1} may have multiplicity higher than one. Appendix C describes the possible cases and provides the generalized eigenvectors when necessary.
6 Characterization of the heteroclinic orbit
The equilibrium M_{0} has a three dimensional unstable manifold and a one dimensional stable manifold while the equilibrium M_{1} has a three dimensional stable manifold and a onedimensional unstable manifold. There is one unstable direction for each equilibrium corresponding to a positive eigenvalue of order {{O}}(\frac{1}{n}). Due to the high dimensionality, it is difficult to read the complete behavior of the flow in phase space. This section aims to develop a picture of the flow on the positive orthant p,q,r,s>0 and to associate the behavior of the system (26) near the singular point \xi=0 with the behavior of the system (S) around M_{0}.
6.1 Behavior near the singular point \xi=0
We begin with the latter point. The following proposition states how (28), (25) are transmitted to the asymptotic behavior of (p,q,r,s) around the equilibrium M_{0} as \eta\to\infty.
Proposition 6.1.
Let \big{(}V,\Sigma,\Theta,\Gamma\big{)}(\xi) be a smooth solution of (26), U(\xi) be defined by (27), and suppose the solution is defined for \xi>0, is smooth, takes values in the positive orthant, and assumes the initial conditions
V(0)=V_{0}\geq 0,\quad\Sigma(0)=\Sigma_{0}>0\,,\quad\Theta(0)=\Theta_{0}>0\,% \quad\Gamma_{0}=\Gamma_{0}>0\quad\mbox{and}\;\;U^{\prime}(0)=0\,. 
Then \big{(}V,\Sigma,\Theta,\Gamma\big{)} and U satisfy at \xi=0 the conditions (29), (33), (34). Morever, the orbit defined by the transformations (37), (39), (42), \chi(\eta)=(p(\eta),q(\eta),r(\eta),s(\eta))\rightarrow M_{0} as \eta\rightarrow\infty. Furthermore, it tends to M_{0} along the direction of the first eigenvector X_{01}, in fact
e^{2\eta}\big{(}\chi(\eta)M_{0}\big{)}\rightarrow\kappa X_{01},\quad\text{% for some constant $\kappa>0$ as $\eta\rightarrow\infty$.}  (46) 
Remark 6.1.
The orbit approaches M_{0} tangent to X_{01} as \eta\rightarrow\infty. Since M_{0} has a threedimensional unstable manifold and \mu_{02}(=1)<\mu_{01}(=2)<\mu_{03}(={{O}}(\frac{1}{n})), the orbits emanating from M_{0} tangent to X_{01} all lie on a two dimensional manifold that at M_{0} is tangent to the plane spanned by X_{01} and X_{03}. This two dimensional submanifold will be referred to as the Strongly unstable manifold of M_{0}.
Proof.
Assuming smoothness and boundedness of \big{(}V,\Sigma,\Theta,\Gamma\big{)} and U in a neighborhood of \xi=0 and the conditions (25), (28) we deduce first (29), (33), (34) by the argument of section 3. The derivatives of \big{(}\Gamma,V,\Theta,\Sigma,U\big{)} at \xi=0 are obtained by differentiating the system (26) repeatedly. Rewrite (26) as
\displaystyle a+\lambda\xi\frac{\Gamma^{\prime}}{\Gamma}  \displaystyle=\frac{U}{\Gamma},  \displaystyle c+\lambda\xi\frac{\Theta^{\prime}}{\Theta}  \displaystyle=\frac{\Sigma\Gamma}{\Theta}\frac{U}{\Gamma},  
\displaystyle(b+\lambda)U+\lambda\xi U^{\prime}(\xi)  \displaystyle=\Sigma^{{}^{\prime\prime}}=\Big{(}\frac{\Sigma\Gamma}{\Theta}% \frac{\Theta}{\Gamma}\Big{)}^{{}^{\prime\prime}},  \displaystyle\frac{\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{{}^{\prime\prime% }}}{\frac{\Sigma\Gamma}{\Theta}}  \displaystyle=(1+m+n)\frac{\Gamma^{{}^{\prime\prime}}}{\Gamma}(1+\alpha)\frac% {\Theta^{{}^{\prime\prime}}}{\Theta}+n\frac{\big{(}\frac{U}{\Gamma}\big{)}^{{}% ^{\prime\prime}}}{\frac{U}{\Gamma}}\,, 
from where after a computation we conclude
\displaystyle\frac{U}{\Gamma}(0)=a=r_{0},  \displaystyle\Big{(}\frac{U}{\Gamma}\Big{)}^{\prime}(0)  \displaystyle=0,  \displaystyle\Big{(}\frac{U}{\Gamma}\Big{)}^{{}^{\prime\prime}}(0)  \displaystyle=\frac{\Gamma(0)}{\Sigma(0)}\frac{2(b+\lambda)r_{0}}{\frac{1s_{% 0}}{\lambda}\frac{n}{r_{0}}\Big{(}\frac{2}{s_{0}}+\frac{r_{0}}{\lambda}\Big{)% }\left(\frac{\frac{1}{\lambda}+2}{\frac{1+\alpha}{\lambda}r_{0}+\frac{2}{s_{0}% }}\right)},  
\displaystyle\frac{\Sigma\Gamma}{\Theta}(0)=\frac{c}{a}=s_{0},  \displaystyle\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{\prime}(0)  \displaystyle=0,  \displaystyle\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{{}^{\prime\prime}}(0)  \displaystyle=\frac{n}{r_{0}}\left(\frac{\frac{1}{\lambda}+2}{\frac{1+\alpha}{% \lambda}r_{0}+\frac{2}{s_{0}}}\right)\Big{(}\frac{U}{\Gamma}\Big{)}^{{}^{% \prime\prime}}(0). 
The Taylor expansions of p(\log\xi), q(\log\xi), r(\log\xi) and s(\log\xi) at \xi=0 are computed using (37), (42), (33) and the relations above,
\displaystyle p(\log\xi)  \displaystyle=\frac{{\bar{\gamma}}}{{\bar{\sigma}}}=\frac{\xi^{a_{1}}\Gamma(% \xi)}{\xi^{d_{1}}\Sigma(\xi)}=\xi^{2}\frac{\Gamma(\xi)}{\Sigma(\xi)}=\xi^{2}% \frac{\Gamma(0)}{\Sigma(0)}+o(\xi^{2})\,,  
\displaystyle q(\log\xi)  \displaystyle=b\frac{{\bar{v}}}{{\bar{\sigma}}}=b\frac{\xi^{b_{1}}V(\xi)}{\xi^% {d_{1}}\Sigma(\xi)}=b\xi\frac{V(\xi)}{\Sigma(\xi)}=b\xi^{2}\frac{U(0)}{\Sigma(% 0)}+o(\xi^{2})=\xi^{2}~{}br_{0}\frac{\Gamma(0)}{\Sigma(0)}+o(\xi^{2})\,,  
\displaystyle r(\log\xi)  \displaystyle=\frac{{\bar{u}}}{{\bar{\gamma}}}=\frac{\xi^{1+b_{1}}U(\xi)}{\xi^% {a_{1}}\Gamma(\xi)}=\frac{U(0)}{\Gamma(0)}+\xi\Big{(}\frac{U}{\Gamma}\Big{)}^{% \prime}(0)+\frac{1}{2}\xi^{2}\Big{(}\frac{U}{\Gamma}\Big{)}^{{}^{\prime\prime}% }(0)+o(\xi^{2})  
\displaystyle=\frac{U}{\Gamma}(0)+\xi^{2}\frac{\Gamma(0)}{\Sigma(0)}\frac{(b+% \lambda)r_{0}}{\frac{1s_{0}}{\lambda}\frac{n}{r_{0}}\Big{(}\frac{2}{s_{0}}+% \frac{r_{0}}{\lambda}\Big{)}\left(\frac{\frac{1}{\lambda}+2}{\frac{1+\alpha}{% \lambda}r_{0}+\frac{2}{s_{0}}}\right)},  
\displaystyle s(\log\xi)  \displaystyle=\frac{{\bar{\sigma}}{\bar{\gamma}}}{{\bar{\theta}}}=\frac{\xi^{a% _{1}+d_{1}}\Sigma(\xi)\Gamma(\xi)}{\xi^{c_{1}}\Theta(\xi)}=\frac{\Sigma\Gamma}% {\Theta}(0)+\xi\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{{}^{\prime}}(0)+% \frac{1}{2}\xi^{2}\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{{}^{\prime\prime}% }(0)+o(\xi^{2})  
\displaystyle=\frac{\Sigma\Gamma}{\Theta}(0)+\xi^{2}n\left(\frac{\big{(}\frac{% 1}{\lambda}+2\big{)}\frac{1}{r_{0}}}{\frac{1+\alpha}{\lambda}r_{0}+\frac{2}{s_% {0}}}\right)\frac{\Gamma(0)}{\Sigma(0)}\frac{(b+\lambda)r_{0}}{\frac{1s_{0}}% {\lambda}\frac{n}{r_{0}}\Big{(}\frac{2}{s_{0}}+\frac{r_{0}}{\lambda}\Big{)}% \left(\frac{\frac{1}{\lambda}+2}{\frac{1+\alpha}{\lambda}r_{0}+\frac{2}{s_{0}}% }\right)}+o(\xi^{2}). 
Therefore, taking note of the eigenvector X_{01} in Appendix C, we conclude
\displaystyle\chi(\log\xi)M_{0}=\big{(}p(\log\xi),q(\log\xi),r(\log\xi),s(% \log\xi)\big{)}M_{0}=\frac{\Gamma(0)}{\Sigma(0)}\xi^{2}X_{01}+o(\xi^{2}), 
which implies (46) since \eta=\log\xi\to\infty as \xi\to 0. ∎
Remark 6.2.
For n small enough, we find that second derivatives have definite signs: \displaystyle\Big{(}\frac{U}{\Gamma}\Big{)}^{{}^{\prime\prime}}(0)<0, \displaystyle\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{{}^{\prime\prime}}(0)<0, and
\displaystyle\frac{\Gamma^{{}^{\prime\prime}}(0)}{\Gamma(0)}  \displaystyle=\frac{1}{2\lambda}\Big{(}\frac{U}{\Gamma}\Big{)}^{{}^{\prime% \prime}}(0)<0,  \displaystyle\frac{\Theta^{{}^{\prime\prime}}(0)}{\Theta(0)}  \displaystyle=\frac{1}{2\lambda}\Big{(}s_{0}\Big{(}\frac{U}{\Gamma}\Big{)}^{{}% ^{\prime\prime}}(0)+r_{0}\Big{(}\frac{\Sigma\Gamma}{\Theta}\Big{)}^{{}^{\prime% \prime}}(0)\Big{)}\Big{)}<0,  (47)  
\displaystyle\frac{U^{{}^{\prime\prime}}(0)}{U(0)}  \displaystyle=\frac{\Gamma^{{}^{\prime\prime}}(0)}{\Gamma(0)}+\frac{\big{(}% \frac{U}{\Gamma}\big{)}^{{}^{\prime\prime}}(0)}{\frac{U}{\Gamma}(0)}<0,  \displaystyle\Sigma^{{}^{\prime\prime}}(0)  \displaystyle=(b+\lambda)U(0)>0. 
6.2 A heteroclinic orbit
The behavior of (26) near the singular point \xi=0 suggests to look for an orbit of (S) emanating from M_{0} in the direction of the Strongly unstable manifold. At the other end point, as \eta\rightarrow\infty, we expect that the variables p, q, r and s equilibrate to a bounded state as \eta\rightarrow\infty. This would imply that the variables {\tilde{\gamma}}, {\tilde{v}}, {\tilde{\theta}}, and {\tilde{\sigma}} grow at most exponentially in \eta as seen from (41) and in turn polynomially in \xi. Therefore, we look for a heteroclinic orbit joining M_{0} with M_{1}, the only two equilibria in the sector r>0 and s>0. That is we expect \chi(\eta)\rightarrow M_{1} as \eta\rightarrow\infty. The heteroclinic orbit will be called \chi(\eta) for the rest of this paper. It will be constructed in the next section.
The two end point behaviors can be interpreted geometrically as well: The end point behavior as \eta\rightarrow\infty specifies a nontrivial submanifold of the unstable manifold of the equilibrium M_{0} from which the orbit emanates. This submanifold will turn out to intersect the stable manifold of M_{1} and the intersection of these two manifold is the heteroclinic orbit we search for.
U_{0}, \Gamma_{0} are selected.
Given the heteroclinic \chi(\eta), this subsection is devoted to adapting the initial data. By data we refer to \big{(}\Gamma_{0},\Theta_{0},\Sigma_{0},U_{0}\big{)}.
6.3 Adapting the initial data
Suppose now that an orbit \chi(\eta) has been constructed that satisfies (S), it emanates from M_{0} in the Strongly unstable manifold, i.e. it satisfies (46), and connects to M_{1}. The orbit corresponds to a set of parameters (\lambda,\alpha,m,n). We proceed to show how the data input \big{(}V_{0},\Gamma_{0},\Theta_{0},\Sigma_{0},U_{0}\big{)} fit under the supposed orbit \chi^{\lambda,\alpha,m,n}(\eta). From (34), (29) we know that V_{0}=0 and
\Theta_{0}=c^{\frac{1}{1+\alpha}}\Gamma_{0}^{\frac{m}{1+\alpha}}U_{0}^{\frac{% 1+n}{1+\alpha}},\quad\Sigma_{0}=c^{\frac{\alpha}{1+\alpha}}\Gamma_{0}^{\frac{m% }{1+\alpha}}U_{0}^{\frac{\alphan}{1+\alpha}}. 
By (35), the rate of growth \lambda determines the ratio \frac{U_{0}}{\Gamma_{0}}
\lambda=\Big{(}\frac{U_{0}}{\Gamma_{0}}\frac{2(1+\alpha)n}{D}\Big{)}\frac{D}% {2(1+\alpha)}.  (48) 
Note finally that the restriction (43) on the growth rate \lambda implies a restriction on the ratio
\displaystyle\frac{2(1+\alpha)n}{D}<\frac{U_{0}}{\Gamma_{0}}  \displaystyle<\frac{2(1+\alpha)n}{D}+\frac{4(1+\alpha)(\alphamn)(1+m)}{D(1+% m+n)^{2}}  (49)  
\displaystyle=\frac{2(1+\alpha)}{1+m+n}\frac{n}{D}\left(\frac{4(1+\alpha)(% \alphamn)}{(1+m+n)^{2}}+1\right). 
It remains to resolve only one degree of freedom. The orbit \chi(\eta) emanating from M_{0} in the direction X_{01} satisfies the asymptotic expansion
\chi(\eta)M_{0}=\kappa_{1}e^{\mu_{01}\eta}X_{01}+\kappa_{3}e^{\mu_{03}\eta}X_% {03}+\text{higherorder terms as $\eta\rightarrow\infty$},  (50) 
for two constants \kappa_{1} and \kappa_{3}. Any reparametrization \chi(\eta\eta_{0}), \eta_{0}\in\mathbb{R}, depicts the same heteroclinic orbit. Define \bar{\chi}(\eta)=\chi(\eta\eta_{0}) and we proceed to select \eta_{0} so as to satisfy the data. Then,
\lim_{\eta\rightarrow\infty}\big{(}\bar{\chi}(\eta)M_{0}\big{)}e^{2\eta}=% \lim_{\eta\rightarrow\infty}\Big{(}\big{(}\chi(\eta\eta_{0})M_{0}\big{)}e^{% 2(\eta\eta_{0})}\Big{)}e^{2\eta_{0}}=e^{2\eta_{0}}\kappa_{1}X_{01}. 
On the other hand, from the proof of Proposition 6.1,
\lim_{\eta\rightarrow\infty}\big{(}\bar{\chi}(\eta)M_{0}\big{)}e^{2\eta}=% \frac{\Gamma_{0}}{\Sigma_{0}}X_{01}\,, 
which dictates we fix the last degree of freedom by selecting
\eta_{0}=\frac{1}{2}\log\left(\frac{\Gamma_{0}}{\Sigma_{0}}\kappa_{1}\right).  (51) 
7 Existence via Geometric theory of singular perturbations
This section is devoted to proving the existence of a heteroclinic orbit \chi(\eta) with limiting behavior as determined in section 6.
Theorem 1.
Let (\alpha,m,n) take values in the range (23). Given \lambda>0 satisfying (43) with n=0, there is n_{0}(\lambda,\alpha,m) such that for n\in[0,n_{0}) and \lambda satisfying (43) the system (S) admits a heteroclinic orbit \chi^{\lambda,\alpha,m,n}(\eta) joining the equilibrium M_{0}^{\lambda,\alpha,m,n} to the equilibrium M_{1}^{\lambda,\alpha,m,n} and satisfying the property
\displaystyle e^{2\eta}\big{(}\chi^{\lambda,\alpha,m,n}(\eta)M_{0}^{\lambda,% \alpha,m,n}\big{)}\rightarrow\kappa X_{01}^{\lambda,\alpha,m,n}\quad\text{as $% \eta\rightarrow\infty$ for some $\kappa\neq 0$}.  (52) 
The heteroclinic orbit \chi^{\lambda,\alpha,m,n}(\eta) is achieved by applying the geometric singular perturbation theory. The presence of the small parameter n>0 in the lefthandside of (S){}_{3} provides a fastslow structure to the system, having r as a fast variable and the rest as slow variables. In the interest of the reader, we present some preliminary information. Experts on geometric singular perturbation theory may wish to proceed directly to Sections 7.2, 7.3.
Recall that (S) accounts for a family of dynamical systems parametrized by (\lambda,\alpha,m,n); the heteroclinic orbit will be achieved respectively for each admissible (\lambda,\alpha,m,n). To simplify notations we suppress the dependence on \lambda, \alpha, and m but retain the dependence on n.
7.1 Invariant manifold theory and geometric singular perturbation theory
We state here some rudiments of the geometric singular perturbation theory from [11, 12]. Fenichel’s persistence theorem is developed in [12, Theorem 9.1]. In the present application the versions [25, Theorem 2.2] and [25, Theorem 3.1] are applied.
Let X be a C^{r} vector field in \mathbb{R}^{d} with r\geq 2 and let \bar{\Lambda}=\Lambda\cup\partial\Lambda be a compact, connected C^{r+1} manifold in \mathbb{R}^{d}. F^{t}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d} denotes the time tmap associated with the vector field X and DF^{t} denotes its differential. \bar{\Lambda} is said to be overflowing invariant under X if for every m\in\bar{\Lambda} and t\leq 0, F^{t}(m)\in\bar{\Lambda} and X is pointing strictly outward on \partial\Lambda. T\mathbb{R}^{d}\bar{\Lambda} denotes the tangent bundle of \mathbb{R}^{d} along \bar{\Lambda} and T\bar{\Lambda} denotes the tangent bundle of \bar{\Lambda}. A subbundle E\subset T\mathbb{R}^{d}\bar{\Lambda} is said to be negatively invariant if E\supset DF^{t}(E) for all t\leq 0.
Let E\subset T\mathbb{R}^{d}\bar{\Lambda} be a subbundle that is negatively invariant and contains T\bar{\Lambda}. Given such E, T\mathbb{R}^{d}\bar{\Lambda} then splits into T\mathbb{R}^{d}\bar{\Lambda}=E\oplus E^{\prime}=T\bar{\Lambda}\oplus N\oplus E% ^{\prime}, where N\subset E is any complement of T\bar{\Lambda} in E and E^{\prime}\subset T\mathbb{R}^{d}\bar{\Lambda} is any complement of E in T\mathbb{R}^{d}\bar{\Lambda}. Next, the subundles are distinguished according to the growth or decay rates of the linearized flow as t\to\infty, following [12]: Let m\in\bar{\Lambda} and v^{0}\in T_{m}\Lambda; w^{0}\in N_{m}; x^{0}\in E^{\prime}_{m}; v^{t}=DF^{t}(m)v^{0}; w^{t}=\pi^{N}DF^{t}(m)w^{0}; x^{t}=\pi^{E^{\prime}}DF^{t}(m)x^{0}, where \pi^{N} and \pi^{E^{\prime}} are bundle projections onto N and E^{\prime} respectively. Define
\displaystyle\nu^{s}(m)  \displaystyle\triangleq\inf\Big{\{}\nu>0\>:\>\frac{1}{x^{t}}=o(\nu^{t})% \quad\text{as $t\rightarrow\infty$}\quad\forall x^{0}\in E^{\prime}_{m}\Big{\}}. 
If \nu^{s}(m)<1, define
\displaystyle\sigma^{s}(m)  \displaystyle\triangleq\inf\Big{\{}\sigma>0\>:\>v^{t}=o(x^{t}^{\sigma})% \quad\text{as $t\rightarrow\infty$}\quad\forall x^{0}\in E^{\prime}_{m},v^{0}% \in T_{m}\Lambda\Big{\}}. 
Next, define
\displaystyle\alpha^{u}(m)  \displaystyle\triangleq\inf\Big{\{}\alpha>0\>:\>w^{t}=o(\alpha^{t})\quad% \text{as $t\rightarrow\infty$}\quad\forall w^{0}\in N_{m}\Big{\}}. 
If \alpha^{u}(m)<1, define
\displaystyle\rho^{u}(m)  \displaystyle\triangleq\inf\Big{\{}\rho>0\>:\>\frac{w^{t}}{v^{t}}=o(\rho% ^{t})\quad\text{as $t\rightarrow\infty$}\quad\forall w^{0}\in N_{m},v^{0}\in T% _{m}\Lambda\Big{\}}. 
If \rho^{u}(m)<1, define
\displaystyle\tau^{u}(m)  \displaystyle\triangleq\inf\Big{\{}\tau>0\>:\>\hat{v}^{t}=o\left(\Big{(}% \frac{v^{t}}{w^{t}}\Big{)}^{\tau}\right)\quad\text{as $t\rightarrow% \infty$}\quad\forall w^{0}\in N_{m},v^{0}\in T_{m}\Lambda,\hat{v}^{0}\in T_{m}% \Lambda\Big{\}}. 
Definition 7.1.
Let \bar{\Lambda}=\Lambda\cup\partial\Lambda be an overflowing invariant manifold such that T\mathbb{R}^{d}\bar{\Lambda} admits a splitting by E as described above. We say an overflowing invariant manifold \bar{\Lambda} satisfies assumptions (A_{r}) and (B_{r^{\prime}}), r^{\prime}\leq r1 if for all m\in\bar{\Lambda} the growth rates hold
\displaystyle\nu^{s}(m)  \displaystyle<1,\quad\sigma^{s}(m)<\frac{1}{r},  (A_{r})  
\displaystyle\alpha^{u}(m)  \displaystyle<1,\quad\rho^{u}(m)<1,\quad\tau^{u}(m)<\frac{1}{r^{\prime}}.  (B_{r^{\prime}}) 
Remark 7.1.
Given the bundle splitting, the conditions (A_{r}) and (B_{r^{\prime}}) suffice to construct the unstable manifold of \bar{\Lambda} as well as the finer foliation structure within it; see [10, Theorem 4] and [11, Theorem 3]. Moreover, for the special case E=T\bar{\Lambda}, Fenichel in [9] proved the persistence of the overflowing manifold \bar{\Lambda} of X_{0} when only (A_{r}) is assumed.
Next, suppose that a C^{r} family of vector fields X_{\epsilon} is given depending on a small parameter for \epsilon\in[\epsilon_{0},\epsilon_{0}] is given. If \bar{\Lambda}_{\epsilon} exists as a C^{r} family of overflowing invariant manifolds for each sufficiently small \epsilon, then the unstable manifold and its foliation structure are persistent in an appropriate sense; see section 16 of [12], and the discussion in [12, p. 90].
Next, we introduce the notion of normal hyperbolicity for manifolds without boundary, [12, p.89] and [9, p.221].
Definition 7.2 (Normally Hyperbolic Invariant Manifold [12]).
Let \Lambda be a compact, invariant under X manifold without boundary. Let E^{s} and E^{u} be subbundles of T\mathbb{R}^{d}\Lambda such that E^{s}+E^{u}=T\mathbb{R}^{d}\Lambda, E^{s}\cap E^{u}=T\Lambda, E^{u} is negatively invariant under X and E^{s} is negatively invariant under X. We say \Lambda is rnormally hyperbolic if \Lambda is an overflowing invariant manifold with a subbundle E^{u} satisfying the rate assumptions (A_{r}) and \Lambda is so with E^{s} under X.
The geometric singular perturbation theory [12, Theorem 9.1] applies the invariant manifold theory to the fastslow structure induced by the dynamical system
\left\{\begin{aligned} \displaystyle\dot{x}&\displaystyle=f(x,y,\epsilon),\\ \displaystyle\epsilon\dot{y}&\displaystyle=g(x,y,\epsilon),\end{aligned}\right% .\quad\text{where $\epsilon\in(\epsilon_{0},\epsilon_{0}),\epsilon_{0}>0$ % small, $x\in\mathbb{R}^{\ell}$, $y\in\mathbb{R}^{k}$, $\ell+k=d$.}  (53) 
We say x is a slow variable and y is a fast variable. We assume that f and g are sufficiently smooth, and the terms slow and fast variable originate from the two limiting asymptotic problems :
\text{(Reduced Problem)}\quad\left\{\begin{aligned} \displaystyle\dot{x}&% \displaystyle=f(x,y,0),\\ \displaystyle 0&\displaystyle=g(x,y,0),\end{aligned}\right.\hskip 25.0pt\text{% (Layer Problem)}\quad\left\{\begin{aligned} \displaystyle x^{\prime}&% \displaystyle=0,\\ \displaystyle y^{\prime}&\displaystyle=g(x,y,0),\quad(\cdot)^{\prime}=\frac{d}% {d(t/\epsilon)}.\end{aligned}\right. 
The zeroset \mathrsfs{S} of g(x,y,0) defines a manifold the orbits of the Reduced problem take values. In general, \mathrsfs{S} is not realized as a graph as it can have many branches. This manifold consists of equilibria of the Layer problem. We consider
\displaystyle\mathrsfs{S}  \displaystyle=\Big{\{}(x,y)\>\Big{}\>g(x,y,0)=0\Big{\}},  
\displaystyle\mathrsfs{S}_{R}  \displaystyle\subset\Big{\{}(x,y)\in\mathrsfs{S}\>\Big{}\>\text{$D_{y}g(x,y,0% )$ has the full rank $k$}\Big{\}}\quad\text{open},  
\displaystyle\mathrsfs{S}_{H}  \displaystyle\subset\Big{\{}(x,y)\in\mathrsfs{S}_{R}\>\Big{}\>\text{all % eigenvalues of $D_{y}g(x,y,0)$ have nontrivial real parts}\Big{\}}\quad\text{% open}. 
On \mathrsfs{S}_{R}, the equation 0=g(x,y,0) is locally solvable for y in terms of x and we speak of the reduced vector field X_{R} on slow variables. (See equation (7.8) in [12].) On a compact subset K\subset\mathrsfs{S}_{H} Fenichel’s persistence Theorem [12, Theorem 9.1] applies.
In the sequel, we need the notion of transversal intersections:
Definition 7.3 (Transversal Intersection).
([25, Definition 3.1]) Let {\mathrsfs{M}}_{1} and {\mathrsfs{M}}_{2} be submanifolds of a manifold {\mathrsfs{M}}. The manifolds {\mathrsfs{M}}_{1} and {\mathrsfs{M}}_{2} intersect transversally at a point m\in{\mathrsfs{M}}_{1}\cap{\mathrsfs{M}}_{2} iff
T_{m}{\mathrsfs{M}}=T_{m}{\mathrsfs{M}}_{1}+T_{m}{\mathrsfs{M}}_{2} 
holds, where T_{m}\mathrsfs{M} denotes the tangent space of the manifold \mathrsfs{M} and similarly for \mathrsfs{M}_{1} and \mathrsfs{M}_{2}.
It is shown in [25] that when heteroclinic orbits are realized as transversal interesections of a stable and an unstable manifold they are persistent under perturbations.
7.2 Singular orbits for the inviscid system with n=0
Let us describe how we apply the perturbation theory to the case of (S). We take two normally hyperbolic manifolds \mathrsfs{N}_{0} and \mathrsfs{N}_{1}, which are simply the equilibrium points M_{0} and M_{1}. The goal of this section is to establish the transversal intersection of the \mathrsfs{N}_{0}^{u}, the unstable manifold of \mathrsfs{N}_{0} in pqrsspace, and \mathrsfs{N}^{s}_{1}, the stable manifold of \mathrsfs{N}_{1}.
A bunch of symbols and series of notations, following [25], are introduced and used for the remaining section. K\subset\mathrsfs{S}_{H} denotes a compact critical manifold to be specified later. The reduced vector field problem is defined in K and X_{R} denotes the reduced vector field. As the reduced vector field is defined in K, K may contain a finer invariant manifold \mathrsfs{N}^{\prime}\hookrightarrow K that is normally hyperbolic to the reduced vector field. As usual \mathrsfs{N^{\prime}}, W^{u}(\mathrsfs{N}^{\prime}), and W^{s}(\mathrsfs{N}^{\prime}) will be respectively the manifold, its unstable, and its stable manifolds embedded in K. In particular, W_{0}^{u} denotes the reduced unstable manifold of \mathrsfs{N}_{0} and W_{1}^{s} denotes the reduced stable manifold of \mathrsfs{N}_{1} embedded in K. Inside of W^{u}(\mathrsfs{N}^{\prime}) are foliations \mathrsfs{F}^{u}_{x} (see [12, Theorem 12.2]), meaning a foliation that passes through x\in W^{u}(\mathrsfs{N}^{\prime}). \mathrsfs{F}^{s}_{x} stands for the analogous foliation of W^{s}(\mathrsfs{N}^{\prime}). Finally, when geometric objects are extended via the invariant manifold theory, those objects for n>0 are denoted by a superscript n, for example K^{n}, \mathrsfs{N}_{0}^{n}, \mathrsfs{N}_{0}^{u,n}, W_{0}^{u,n}, M_{0}^{n}, \cdots.
Next, we identify the Reduced problem,
\displaystyle\dot{p}  \displaystyle=p\Big{(}\frac{1}{\lambda}({r}a)+2\lambda p{r}q\Big{)},  (R)  
\displaystyle\dot{q}  \displaystyle=q\Big{(}1\lambda p{r}q\Big{)}+bp{r},  
\displaystyle 0  \displaystyle=r\Big{(}\frac{\alpham}{\lambda(1+\alpha)}(ra)+\lambda pr+q+% \frac{\alpha}{\lambda}r\big{(}s\frac{1+m}{1+\alpha}\big{)}\Big{)},  
\displaystyle\dot{s}  \displaystyle=s\Big{(}\frac{\alpham}{\lambda(1+\alpha)}({r}a)+\lambda p{r}+q% \frac{1}{\lambda}{r}\big{(}s\frac{1+m}{1+\alpha}\big{)}\Big{)}, 
and the Layer problem for the system (S),
\displaystyle{p}^{\prime}  \displaystyle=0,\quad{q}^{\prime}=0,\quad{r}^{\prime}=r\Big{(}\frac{\alpham}{% \lambda(1+\alpha)}(ra)+\lambda pr+q+\frac{\alpha}{\lambda}r\big{(}s\frac{1+m% }{1+\alpha}\big{)}\Big{)},\quad{s}^{\prime}=0.  (54) 
Here, (\cdot)^{\prime}=\frac{d}{d\tilde{\eta}}:=\frac{d}{d(\eta/n)} denotes differentiation with respect to the fast independent variable \tilde{\eta}. The zeroset of the function
g(p,q,r,s)\triangleq r\Big{(}\frac{\alpham}{\lambda(1+\alpha)}(ra)+\lambda pr% +q+\frac{\alpha}{\lambda}r\big{(}s\frac{1+m}{1+\alpha}\big{)}\Big{)}  (55) 
consists of the equilibria of (54).
7.2.1 Choice of the critical manifold K
The algebraic equation g(p,q,r,s)=0 specifies three dimensional hypersurfaces. Away from the r\equiv 0 plane, one may obtain the hypersurface as a graph of the function
r=\hat{r}(p,q,s)=\frac{\frac{\alpham}{\lambda(1+\alpha)}aq}{\frac{\alpham}{% \lambda(1+\alpha)}+\lambda p+\frac{\alpha}{\lambda}\big{(}s\frac{1+m}{1+% \alpha}\big{)}}\,, 
or implicitly in the form
\frac{\alpham}{\lambda(1+\alpha)}(\hat{r}a)+\lambda p\hat{r}+q+\frac{\alpha}% {\lambda}\hat{r}\big{(}s\frac{1+m}{1+\alpha}\big{)}=0.  (56) 
The level sets of \hat{r} can be used in the pqsspace in order to visualize the hypersurface; they are affine due to (56). Fig. 3 illustrates a few marked affine surfaces of \hat{r}(p,q,s)=R, in the range 0\leq R\leq a. When R=a(=r_{0}), it passes through \left(0,0,\tfrac{1+m}{1+\alpha}\right), which is the equilibrium M_{0}. As R decreases the affine level sets sweep the positive p,q sector. The surface crosses the other equilibrium M_{1} when R=r_{1}. Then R continues to decrease until it touches the r\equiv 0 plane.
The critical manifold K is selected taking account of the properties of \hat{r}(p,q,r). The domain of \hat{r} is a trapezoid in pqsspace
\displaystyle D  \displaystyle\triangleq\left\{\>(p,q,s)\>\Big{}\>p\geq\epsilon,~{}~{}q\leq 2% ,~{}~{}\lefts\frac{1+m}{1+\alpha}\right\leq\frac{1}{2}\min\left\{\frac{% \alpham}{\alpha(1+\alpha)},\frac{1+m}{(1+\alpha)}\right\},\right.  
\displaystyle\left.\hat{r}(p,q,s)\geq\frac{1}{2}\min\{1,r_{1}\}\right\}. 
where \epsilon is a positive parameter selected sufficiently small. K is then defined by setting K\triangleq\big{(}D,\hat{r}(D)\big{)}.
Note that K is chosen so that (i) M_{0} and M_{1} are on K; (ii) s and r=\hat{r}(p,q,s) have positive lower bound on K. See the trapezoid D in Fig. 4.
Next, we verify that K\subset\mathrsfs{S}_{H}.
Proposition 7.1.
We have that K\subset\mathrsfs{S}_{H}, i.e., the partial jacobian \frac{\partial g}{\partial r}(p,q,r,s)_{r=\hat{r}(p,q,s)}>0 for all (p,q,r,s)\in K.
Proof.
\displaystyle\left.\frac{\partial g}{\partial r}\right_{K}  \displaystyle=\Big{(}\frac{\alpham}{\lambda(1+\alpha)}(\hat{r}a)+\lambda p% \hat{r}+q+\frac{\alpha}{\lambda}\hat{r}\big{(}s\frac{1+m}{1+\alpha}\big{)}% \Big{)}+\hat{r}\Big{(}\frac{\alpham}{\lambda(1+\alpha)}+\lambda p+\frac{% \alpha}{\lambda}\big{(}s\frac{1+m}{1+\alpha}\big{)}\Big{)}  
\displaystyle=\hat{r}\Big{(}\frac{\alpham}{\lambda(1+\alpha)}+\lambda p+\frac% {\alpha}{\lambda}\big{(}s\frac{1+m}{1+\alpha}\big{)}\Big{)}\geq\frac{1}{2}% \min\{1,r_{1}\}\Big{(}\frac{\alpham}{2\lambda(1+\alpha)}\lambda\epsilon\Big{% )}. 
It suffices to take \epsilon<\frac{\alpham}{4\lambda^{2}(1+\alpha)}, independently of n. ∎
7.2.2 Nested invariant manifold structures in K
The flow (R), strictly restricted on K, is further analyzed. The three dimensional flow
\displaystyle\dot{p}  \displaystyle=p\Big{(}\frac{D}{\lambda(1+\alpha)}(\hat{r}a_{0})+\frac{\alpha}% {\lambda}\hat{r}\big{(}s\frac{1+m}{1+\alpha}\big{)}\Big{)},  (\text{R}^{\prime})  
\displaystyle\dot{q}  \displaystyle=q\Big{(}1\lambda p\hat{r}q\Big{)}+bp\hat{r},  
\displaystyle\dot{s}  \displaystyle=\frac{1+\alpha}{\lambda}\hat{r}s\Big{(}s\frac{1+m}{1+\alpha}% \Big{)}. 
augmented in the rdirection by r=\hat{r}(p,q,s) is the flow of the Reduced system (R).
It is necessary to pinpoint a few finer invariant structures of the reduced flow (\text{R}^{\prime}), on which the invariant manifold theory can equally well be applied. This will be crucial ingredient of our arguments. To summarize, in the three dimensional reduced space K, we consider the embeddings
M_{0},M_{1}\quad\hookrightarrow\quad\text{$p\equiv 0$ line on $s\equiv\tfrac{1% +m}{1+\alpha}$ plane}\quad\hookrightarrow\quad\text{$s\equiv\tfrac{1+m}{1+% \alpha}$ plane}\quad\hookrightarrow\quad K\,, 
which consist of manifolds all of which are invariant under the flow (\text{R}^{\prime}). Indeed, on the plane s\equiv\tfrac{1+m}{1+\alpha}, (\text{R}^{\prime}) decouples,
\displaystyle\dot{p}  \displaystyle=p\Big{(}\frac{D}{\lambda(1+\alpha)}(\hat{r}a_{0})\Big{)},  (57)  
\displaystyle\dot{q}  \displaystyle=q\Big{(}1\lambda p\hat{r}q\Big{)}+bp\hat{r}, 
where \hat{r}=\hat{r}\big{(}p,q,\frac{1+m}{1+\alpha}\big{)}. Restricting to p\equiv 0 we obtain yet another invariant line and importantly this line contains the equilibrium points M_{0} and M_{1}.
Fig. 5 illustrates the rest of the program. The justification of the following descriptions will be the subject of the next section. M_{1} is a stable node; the three dimensional volume surrounding M_{1} in Fig. 5 depicts its stable manifold. M_{0} is a saddle; M_{0} has two unstable dimensions in s\equiv\tfrac{1+m}{1+\alpha}, and has one stable dimension in its oblique direction. The vector aligned to qaxis and the one to the green orbit are two eigenvectors for the unstable dimensions. This explains how our heteroclinic orbit (the green one) appears in K.
Not all of the manifolds appearing are normally hyperbolic to the reduced vector field (\text{R}^{\prime}): M_{0} and M_{1} are hyperbolic equilibrium points; \hat{\mathrsfs{N}}^{0}, a segment of p\equiv 0 line (the blue portion in Fig. 5) will be identified as an overflowing manifold satisfying the rate assumptions (A_{r}) and (B_{r^{\prime}}). However, the plane s\equiv\tfrac{1+m}{1+\alpha}, in general, does not satisfy the necessary rate assumptions. See Remark 7.2
7.2.3 Analysis of the flow of (\text{R}^{\prime}).
In the phase space K, the flow of (\text{R}^{\prime}) can be completely analyzed. We visualize the overall flow by first analyzing the flow when restricted to the invariant plane s\equiv\tfrac{1+m}{1+\alpha}, and then noting that off the invariant plane the flow amounts to a stable relaxation process towards the invariant plane at s\equiv\tfrac{1+m}{1+\alpha}, see Fig. 6. The flow in the plane s\equiv\tfrac{1+m}{1+\alpha} is characterized by using planar dynamical systems theory.
Remark 7.2.
The property that the invariant plane s\equiv\tfrac{1+m}{1+\alpha} has one stable direction in K may lead one to believe that the plane is normally hyperbolic. Indeed, the plane s\equiv\tfrac{1+m}{1+\alpha} admits a splitting T\mathbb{R}^{3}=T\bar{\Lambda}\oplus E^{s} with the normal direction decaying. However, the notion of normal hyperbolicity requires stronger properties than merely admitting a splitting: Note that in (A_{r}) and (B_{r^{\prime}}) upon a given splitting we demand \sigma^{s}<1/r and \rho^{u}<1/r^{\prime}. The latter for example leads to demanding that the negative eigenvalue \mu_{14} of M_{1} is strictly less than all other negative eigenvalues, which is not always the case. As a consequence, the persistence Theorem does not apply and we do not assert its persistence under perturbations.
We turn to the reduced linear stability of M_{0} and M_{1} in pqsspace. To summarize, the expressions in Section 5 also hold for n=0, with the exception that the third eigenvalue \mu_{03} and the third eigenvector X_{03} of M_{0} (\mu_{13} and X_{13} respectively of M_{1}) are no longer used. Indeed when n=0, the flow is restricted on the three dimensional set K.
The following Lemma utilizes planar dynamical systems theory to characterizes the flow on a triangle T in the pqplane.
Lemma 7.1.
Let T be the closed triangle on s\equiv\frac{1+m}{1+\alpha} enclosed by p=0, q=0, and the level set \hat{r}(p,q,\frac{1+m}{1+\alpha})=\frac{1}{2}\min\{1,r_{1}\} that is intersected by s\equiv\frac{1+m}{1+\alpha}. Then T\setminus M_{0}\subset W^{s}(M_{1}).
Proof.
T is a two dimensional compact positively invariant set: (1) on p=0, \nu=(1,0) and X_{R}\cdot\nobreak\nu=\dot{p}=0, where X_{R} stands for the reduced vector field of (\text{R}^{\prime}); (2) on q=0, \nu=(0,1) and X_{R}\cdot\nobreak\nu=\dot{q}=bp\hat{r}\geq 0 (b in (21) is always positive); lastly on the hypotenuse, let \underline{r}=\frac{1}{2}\min\{1,r_{1}\}. The inward normal vector is \nu=(\lambda\underline{r},1). We compute
\displaystyle X_{R}\cdot\nu=\lambda\underline{r}\dot{p}\dot{q}  \displaystyle=\lambda\underline{r}p\Big{(}1\lambda\underline{r}pq+\frac{1}{% \lambda}(\underline{r}a)+1\Big{)}q(1\lambda\underline{r}pq\big{)}b% \underline{r}p  
\displaystyle=(1\lambda\underline{r}pq)(\lambda\underline{r}pq)\underline% {r}p\Big{(}(\underline{r}a)+\lambda+b\Big{)}  
\displaystyle=\left(\frac{\alpham}{\lambda(1+\alpha)}\right)^{2}(\underline{r% }r_{0})(\underline{r}r_{1})+\underline{r}p(1\underline{r})\geq\delta>0.  (58) 
Let \Omega be an \omegalimit set of the orbit from x_{0}\in T\setminus M_{0}. It is nonempty because T is compact. It cannot contain M_{0}, because for the flow restricted in T, M_{0} does not have a stable manifold. It cannot contain a periodic orbit; if it did then there would be a fixed point in the interior of T and this is not the case. It cannot contain a separatrix cycle because T has only two fixed points M_{0} and M_{1} and again M_{0} does not have a stable manifold. By PoincaréBendixson Theorem, the \omegalimit set is M_{1}. ∎
Now we are able to state: We call \mathrsfs{F}^{u}_{M_{0}}\subset W^{u}_{0} the strongly unstable manifold of M_{0} satisfying (52) (the green line in Fig. 5) that is characterizable by the Unstable manifold theorem for the hyperbolic fixed point. That \mathrsfs{F}^{u}_{M_{0}} ends up arriving at M_{1} follows by Lemma 7.1, and this gives the proof for n=0 of Theorem 1. The following proposition shows that the one dimensional manifold \mathrsfs{F}^{u}_{M_{0}}\subset W^{u}_{0} intersects the three dimensional manifold W_{1}^{s}(=W^{s}(M_{1})) transversally (see Fig. 5).
Proposition 7.2.
Let \mathrsfs{N}_{0}=M_{0}, \mathrsfs{N}_{1}=M_{1}, \mathrsfs{F}^{u}_{M_{0}}\subset W^{u}_{0} the strongly unstable manifold of M_{0} satisfying (52), W^{s}_{1}=\Phi_{t_{0}}(W^{s}_{loc}(M_{1})), the time t_{0} image of the local stable manifold of M_{1} for large enough t_{0}<\infty. Then \mathrsfs{F}^{u}_{M_{0}} intersects W^{s}_{1} transversally in pqsspace.
Proof of Proposition 7.2.
For large enough t_{0}<\infty, by Lemma 7.1 the orbit point x\in\mathrsfs{F}^{u}_{M_{0}} must be attained in W^{s}_{1} as an interior point. Therefore the tangent space {T}_{x}W^{s}_{1} is the whole of {T}_{x}\mathbb{R}^{3}. Then the intersection with \mathrsfs{F}^{u}_{M_{0}} is trivially transversal. ∎
7.3 Persistence for n>0
Having set forth the critical manifold K in \mathrsfs{S}_{H} and the reduced vector field X_{R} on K, the theorem of Fenichel holds in K; the family K^{n} of slow manifolds persistently exist provided n is sufficiently small. Now we show the finer hyperbolic structure of \mathrsfs{F}^{u}_{M_{0}}\hookrightarrow K.
Lemma 7.2.
Let \mathrsfs{N}_{0}=M_{0}, \mathrsfs{F}^{u}_{M_{0}}\subset W_{0}^{u} the strongly unstable manifold of M_{0} satisfying (52). Then, for sufficiently small n, \mathrsfs{F}^{u}_{M_{0}} perturbs in a C^{r1} manner to \mathrsfs{F}^{u,n}_{M_{0}^{n}} the strongly unstable manifold of M_{0}^{n} satisfying (52).
Proof.
In the pqsspace we select the line segment \hat{\mathrsfs{N}}^{0} that is the transversal intersection of the invariant plane \Big{\{}(p,q,s)~{}~{}p=0\text{ and }q\in[\frac{1}{2},\frac{1}{2}]\Big{\}} and the unstable manifold W_{0}^{u}. \hat{\mathrsfs{N}}^{n} is defined as the intersection of the same plane with W_{0}^{u,n}. We see that \hat{\mathrsfs{N}}^{0}=\Big{\{}(p,q,s)~{}~{}p=0,~{}q\in[\frac{1}{2},\frac{1}% {2}]\text{ and }s=\frac{1+m}{1+\alpha}\Big{\}}.
\hat{\mathrsfs{N}}^{0} is the one dimensional orbit in W_{0}^{u} that is not strongly unstable. We claim that \hat{\mathrsfs{N}}^{0} is an overflowing invariant manifold as in Definition 7.1 of the reduced problem. More precisely, it satisfies (A_{r}) and (B_{r^{\prime}}) with r^{\prime}=r1 and E the tangent pqplane.
From (57), \dot{q}=q(1q) on qaxis, it is clear that \hat{\mathrsfs{N}}^{0} is overflowing invariant. Let E be pqplane along \hat{\mathrsfs{N}}^{0} and E^{\prime} be the lines parallel to the saxis. Then, T\mathbb{R}^{3}\hat{\mathrsfs{N}}^{0} splits into three one dimensional bundles T\hat{\mathrsfs{N}}^{0}\oplus N\oplus E^{\prime} with N complementary to T\hat{\mathrsfs{N}}^{0} in E such that N_{M_{0}} is parallel to X_{01}. The asymptotic rates are determined at M_{0} by the eigenvalues of M_{0}. At M_{0}, E^{\prime}_{M_{0}} is the stable subspace with eigenvalue \mu_{04} and N_{M_{0}} and T_{M_{0}}\hat{\mathrsfs{N}}^{0} are the unstable ones with \mu_{01}=2 and \mu_{02}=1 respectively. From these, we compute
\nu^{s}=e^{\mu_{04}},\quad\sigma^{s}=0,\quad\alpha^{u}=e^{2},\quad\rho^{u}=e% ^{1},\quad\tau^{u}=0. 
Therefore, for the given family of overflowing manifolds \hat{\mathrsfs{N}}^{n}, the strongly unstable manifold and its foliations \mathrsfs{F}^{u}(x,n):=\mathrsfs{F}^{u,n}_{x} exist as a C^{r1} family in both arguments x and n. In turn, the foliation \mathrsfs{F}^{u}(M_{0}^{n},n) that passes through M_{0}^{n} is a C^{r1} map in n. ∎
Persistence of the stable manifold W_{1}^{s}(=W^{s}(M_{1})) is a consequence of the classical stable manifold theorem. Theorem 1 follows in the same way as in [25, Theorem 3.1] by the transversal intersection.
Proof of Theorem 1.
By the theorem of Fenichel, for given (\lambda,\alpha,m,0) satisfying (23) and (48), n_{0} can be taken sufficiently small so that if n\in[0,n_{0}) then (\lambda,\alpha,m,n) satisfies (23) and (48) and the system (S) admits a transversal heteroclinic orbit joining equilibrium M_{0}^{n} to equilibrium M_{1}^{n}: \mathrsfs{F}^{u}_{M_{0}} perturbs to \mathrsfs{F}^{u,n}_{M_{0}^{n}} by Lemma 7.2 and W_{1}^{s} perturbs to W_{1}^{s,n} and the transversal intersection is stable under the perturbation. ∎
8 Emergence of localization
By transforming back using (24), (37), (39), and (42), we recover the profile \big{(}\Gamma(\xi),V(\xi),\Theta(\xi),\Sigma(\xi)\big{)} and U(\xi) by (27) and the associated solution. We replace t\rightarrow t+1 to obtain the final expression:
\displaystyle\gamma(t,x)  \displaystyle=(t+1)^{a}\Gamma((t+1)^{\lambda}x),  \displaystyle v(t,x)  \displaystyle=(t+1)^{b}V((t+1)^{\lambda}x),  \displaystyle\theta(t,x)  \displaystyle=(t+1)^{c}\Theta((t+1)^{\lambda}x),  
\displaystyle\sigma(t,x)  \displaystyle=(t+1)^{d}\Sigma((t+1)^{\lambda}x),  \displaystyle u(t,x)  \displaystyle=(t+1)^{b+\lambda}U((t+1)^{\lambda}x)\,. 
We interpret \big{(}\Gamma(\xi),V(\xi),\Theta(\xi),\Sigma(\xi)\big{)}=\big{(}\gamma(0,x),v(% 0,x),\theta(0,x),\sigma(0,x)\big{)}_{x=\xi} as the initial state. For given material parameters (\alpha,m,n), there are two available degrees of freedom giving rise to a twoparameters family of solutions. As noted in Section 6.3, the choices of U_{0} and \Gamma_{0} determine the selfsimilar profile while the remaining boundary values (\Theta_{0}, \Sigma_{0}) and the rate \lambda are induced by them. The range of U_{0} and \Gamma_{0} is such that
\frac{2(1+\alpha)n}{D}<\frac{U_{0}}{\Gamma_{0}}<\frac{2(1+\alpha)n}{D}+\frac% {4(1+\alpha)(\alphamn)(1+m)}{D(1+m+n)^{2}}. 
The localizing rate \lambda satisfies (48) and takes values 0<\lambda<\frac{2(\alphamn)}{1+m+n}\left(\frac{1+m}{1+m+n}\right). In the sequel, we establish properties of the profiles and the emergence of localization, in the sense of definition (5), (6).
8.1 Properties of the selfsimilar profiles
We first list some information on the behavior of the profiles near \xi=0 and as \xi\to\infty. The latter determines the behavior of the induced solutions off the localization zone.
Proposition 8.1.
Let \big{(}\Gamma(\xi),V(\xi),\Theta(\xi),\Sigma(\xi)\big{)} be the selfsimilar profiles defined by transformations of (37), (39), (42) from the heteroclinic orbit \chi(\eta)=\big{(}p(\eta),q(\eta),r(\eta),s(\eta)\big{)} constructed in Theorem 1 in the range of parameters \Gamma(0)=\Gamma_{0} and U(0)=U_{0} depicted by (49). U(\xi) is defined by (27). Then,

The selfsimilar profile achieves the boundary condition at \xi=0,
{V}(0)=\Gamma_{\xi}(0)=\Theta_{\xi}(0)=\Sigma_{\xi}(0)={U}_{\xi}(0)=0,\quad% \Gamma(0)=\Gamma_{0},\quad U(0)=U_{0}. 
Its asymptotic behavior as \xi\rightarrow 0 is given by
\displaystyle\Gamma(\xi)\Gamma_{0} \displaystyle=\Gamma^{{}^{\prime\prime}}(0)\frac{\xi^{2}}{2}+o(\xi^{2}), \displaystyle\Gamma^{{}^{\prime\prime}}(0) \displaystyle<0, (59) \displaystyle\Theta(\xi)c^{\frac{1}{1+\alpha}}\Gamma_{0}^{\frac{m}{1+\alpha}% }U_{0}^{\frac{1+n}{1+\alpha}} \displaystyle=\Theta^{{}^{\prime\prime}}(0)\frac{\xi^{2}}{2}+o(\xi^{2}), \displaystyle\Theta^{{}^{\prime\prime}}(0) \displaystyle<0, \displaystyle\Sigma(\xi)c^{\frac{\alpha}{1+\alpha}}\Gamma_{0}^{\frac{m}{1+% \alpha}}U_{0}^{\frac{\alphan}{1+\alpha}} \displaystyle=\Sigma^{{}^{\prime\prime}}(0)\frac{\xi^{2}}{2}+o(\xi^{2}), \displaystyle\Sigma^{{}^{\prime\prime}}(0) \displaystyle>0, \displaystyle U(\xi)U_{0} \displaystyle=U^{{}^{\prime\prime}}(0)\frac{\xi^{2}}{2}+o(\xi^{2}), \displaystyle U^{{}^{\prime\prime}}(0) \displaystyle<0, \displaystyle V(\xi)U_{0}\xi \displaystyle=U^{{}^{\prime\prime}}(0)\frac{\xi^{3}}{6}+o(\xi^{3}), \displaystyle U^{{}^{\prime\prime}}(0) \displaystyle<0. 
Its asymptotic behavior as \xi\rightarrow\infty is given by
if \mu_{11}\neq1, or \mu_{11}=1 but b=\lambda,\displaystyle\Gamma(\xi) \displaystyle={{O}}\big{(}\xi^{\frac{1+\alpha}{\alphamn}}), \displaystyle V(\xi) \displaystyle={{O}}\big{(}1), \displaystyle\Theta(\xi) \displaystyle={{O}}\big{(}\xi^{\frac{1+m+n}{\alphamn}}), (60) \displaystyle\Sigma(\xi) \displaystyle={{O}}\big{(}\xi), \displaystyle U(\xi) \displaystyle={{O}}\big{(}\xi^{\frac{1+\alpha}{\alphamn}}) otherwise
\displaystyle\Gamma(\xi) \displaystyle={{O}}\big{(}\xi^{\frac{1+\alpha}{\alphamn}}\big{(}\log\xi\big% {)}^{\frac{1+\alpha}{D}}\big{)}, \displaystyle V(\xi) \displaystyle={{O}}\big{(}\big{(}\log\xi\big{)}^{\frac{\alphamn}{D}}\big{)}, (61) \displaystyle\Theta(\xi) \displaystyle={{O}}\big{(}\xi^{\frac{1+m+n}{\alphamn}}\big{(}\log\xi\big{)}% ^{\frac{1+m+n}{D}}\big{)}, \displaystyle\Sigma(\xi) \displaystyle={{O}}\big{(}\xi\big{(}\log\xi\big{)}^{\frac{\alphamn}{D}}\big% {)}, \displaystyle U(\xi) \displaystyle={{O}}\big{(}\xi^{\frac{1+\alpha}{\alphamn}}\big{(}\log\xi\big% {)}^{\frac{1+\alpha}{D}}\big{)}
Proof.
The proof of the Proposition 6.1 and Remark 6.2 contains (i) and (ii) and thus we are left to prove (iii). In a similar fashion to (50), any orbit \psi(\eta) in the local stable manifold of W^{s}(M_{1}) is characterized by a triple (\kappa_{1}^{\prime},\kappa_{2}^{\prime},\kappa_{3}^{\prime}) in association with the asymptotic expansion
\displaystyle\psi(\eta)M_{1}  (62)  
\displaystyle=\begin{cases}\kappa_{1}^{\prime}e^{\mu_{11}\eta}X_{11}+\kappa_{2% }^{\prime}e^{\mu_{12}\eta}X_{12}+\kappa_{4}^{\prime}e^{\mu_{14}\eta}X_{14}+% \text{high order terms}&\text{if $\mu_{11}\neq1$, or $\mu_{11}=1$ but $b=% \lambda$,}\\ \kappa_{1}^{\prime}\eta e^{\mu_{11}\eta}X_{11}^{\prime}+\kappa_{2}^{\prime}e^{% \mu_{12}\eta}X_{12}+\kappa_{4}^{\prime}e^{\mu_{14}\eta}X_{14}+\text{high order% terms}&\text{if $\mu_{11}=\mu_{12}=1$ and $b\neq\lambda$}\end{cases} 
as \eta\rightarrow\infty. The second formula reflects the presence of a generalized eigenvector.
Now, we have q\rightarrow 1, r\rightarrow r_{1}, s\rightarrow s_{1} but p\rightarrow 0 and the leading order of p is to be found. We can determine the coefficient of X_{11}, above, because the pcomponent of the vectors X_{12} and X_{14} is 0. Since the plane p\equiv 0 is an invariant plane for a nonlinear flow, triplets of the form (0,\kappa_{2}^{\prime},\kappa_{3}^{\prime}) spans this invariant plane. Because our heteroclinic orbit \chi(\eta) ventures out from the plane p\equiv 0, \kappa_{1}^{\prime} for the expansion of \chi(\eta) cannot be 0. This implies that the leading order of p(\log\xi) is
p(\log\xi)=\begin{cases}{{O}}(\xi^{\mu_{11}})&\text{if $\mu_{11}\neq1$ or $% \mu_{11}=1$ but $b=\lambda$,}\\ {{O}}(\xi^{\mu_{11}}\log\xi)&\text{otherwise}\end{cases} 
as \xi\rightarrow\infty.
Asymptotics (60) and (61) are the straightforward calculations obtained from the reconstruction formulas
\displaystyle{\tilde{\gamma}}  \displaystyle=p^{\frac{1+\alpha}{D}}r^{\frac{n}{D}}s^{\frac{\alpha}{D}},  \displaystyle{\tilde{v}}  \displaystyle=\frac{1}{b}p^{\frac{\alphamn}{D}}qr^{\frac{n}{D}}s^{\frac{% \alpha}{D}},  \displaystyle{\tilde{\theta}}  \displaystyle=p^{\frac{1+m+n}{D}}r^{\frac{2n}{D}}s^{\frac{1mn}{D}},  
\displaystyle{\tilde{\sigma}}  \displaystyle=p^{\frac{\alphamn}{D}}r^{\frac{n}{D}}s^{\frac{\alpha}{D}},  \displaystyle{\tilde{u}}  \displaystyle=p^{\frac{1+\alpha}{D}}r^{\frac{n}{D}+1}s^{\frac{\alpha}{D}}, 
8.2 Emergence of localization
As time proceeds, the initial nonuniformity evolves into localization. This section is devoted to describing the behavior of the various fields as time advances. We only present the generic case \frac{1+m+n}{\alphamn}\neq1; in the nongeneric case we would have to add a logarithmic correction according to Proposition 8.1.

Strain : The strain keeps increasing as time proceeds. The growth at the origin is faster than the growth rate at all other points:
\displaystyle\gamma(t,0) \displaystyle=(1+t)^{\frac{2+2\alphan}{D}+\frac{2+2\alpha}{D}\lambda}\Gamma(0), \displaystyle\gamma(t,x) \displaystyle\sim t^{\frac{2+2\alphan}{D}\frac{(1+\alpha)(1+m+n)}{D(\alpham% n)}\lambda}x^{\frac{1+\alpha}{\alphamn}},\quad\text{as $t\rightarrow% \infty$, $x\neq 0$.} Recall that the condition \frac{2+2\alphan}{D}\frac{(1+\alpha)(1+m+n)}{D(\alphamn)}\lambda>0 was the ground for imposing (43), placed to guarantee that the plastic strain is growing even outside the localization zone. On the other hand, the difference between the rate of growth of \gamma at x=0 and the rate at x\neq 0 is easily computed as \frac{1+\alpha}{\alphamn}\lambda>0, which indicates localization of the profile of \gamma around x=0.

Temperature : For the temperature , the growth at the origin is again faster than other points,
\displaystyle\theta(t,0) \displaystyle=(1+t)^{\frac{2(1+m)}{D}+\frac{2(1+m+n)}{D}\lambda}\Theta(0), \displaystyle\theta(t,x) \displaystyle\sim t^{\frac{2(1+m)}{D}\frac{(1+m+n)^{2}}{D(\alphamn)}\lambda% }x^{\frac{1+m+n}{\alphamn}},\quad\text{as $t\rightarrow\infty$, $x\neq 0$.} Again, the positivity of the growth rate \frac{2(1+m)}{D}\frac{(1+m+n)^{2}}{D(\alphamn)}\lambda is a consequence of (43).

Strain rate : The growth rates of the strainrate is by definition less by one to those of the strain, again illustrating localization.
\displaystyle u(t,0) \displaystyle=(1+t)^{\frac{1+m}{D}+\frac{2+2\alpha}{D}\lambda}U(0), \displaystyle u(t,x) \displaystyle\sim t^{\frac{1+m}{D}\frac{(1+\alpha)(1+m+n)}{D(\alphamn)}% \lambda}x^{\frac{1+\alpha}{\alphamn}},\quad\text{as $t\rightarrow\infty$,% $x\neq 0$.} 
Stress : The stress decays with time at all points, but the decay at x=0 is much faster than the decay in other places, indicating stresscollapse in the interior of the band:
\displaystyle\sigma(t,0) \displaystyle=(1+t)^{\frac{2\alpha+2m+n}{D}+\frac{2\alpha+2m+2n}{D}\lambda}% \Sigma(0), \displaystyle\sigma(t,x) \displaystyle\sim t^{\frac{2\alpha+2m+n}{D}+\frac{1+m+n}{D}\lambda}x,\quad% \text{as $t\rightarrow\infty$, $x\neq 0$,} The difference of the two rates is \big{(}\frac{1+m+n}{D}\frac{2\alpha+2m+n}{D}\big{)}\lambda=\lambda.

Velocity : The velocity is an odd function of x. At fixed t, v(t,x) is an increasing function of x ranging from v_{\infty}(t) to v_{\infty}(t), where v_{\infty}(t)\triangleq\lim_{x\rightarrow\infty}v(t,x). The velocity field is contrasted with the linear field of uniform shear motion. The selfsimilar scaling \xi=(1+t)^{\lambda}x implies that most of the transition takes place around the origin leading eventually to stepfunction behavior as time goes to infinity. The asymptotic velocity is
v_{\infty}(t)=(1+t)^{b}V_{\infty}=(1+t)^{\frac{1+m}{D}+\frac{1+m+n}{D}\lambda}% V_{\infty},\quad V_{\infty}\triangleq\lim_{\xi\rightarrow\infty}V(\xi)<\infty. Note that the far field loading condition is different from the linear profile of uniform shearing. This deviation is a consequence of our simplifying assumption of selfsimilarity.
9 Numerical computation of the heteroclinic orbit
In this section we present in detail the process we followed to capture numerically the heteroclinic orbit connecting M_{0} and M_{1}. This is a challenging computational task since both M_{0} and M_{1} are saddle points and the heteroclinic orbit connecting them is the intersection of two 3dimensional manifolds W^{u}(M_{0}) and W^{s}(M_{1}) in \mathbb{R}^{4}. Here, we use the software package AUTO, [6], [7], [8] to compute the heteroclinic orbit connecting M_{0} and M_{1}. One of the main capabilities of AUTO is that it can perform limited bifurcation analysis for parametric systems of ordinary differential equations of the form : \displaystyle u^{\prime}(t)=f(u(t),\chi) where f(\cdot,\cdot),\ u(\cdot)\in\mathbb{R}^{d} and \chi could be one or a set of free parameters.
A direct application of AUTO for solving system (S) and computing the desired heteroclinic orbit will fail. A more careful approach has to be considered, starting from some well prepared data and continuing by exploiting the continuation capabilities of AUTO. Indeed, we start with an exact solution of system (S) available for a specific value of the parameter \alpha and variable p, followed by a projection and two continuation steps escaping from these particular choices for \alpha and p allowing us to compute the heteroclinic orbit. We proceed by describing these four steps in detail.
9.1 Continuation by AUTO
Step 1. (Exact solution) The system (S) admits an explicit solution for certain values of the parameters \alpha,n and the variable p. First for \alpha=0 the equations for p,q,r in (S) decouple from the equation of s. This reduced system carries only the three dimensional unstable manifold of M_{0} characterizing the heteroclinic orbit as a nodesaddle connection. Hence, in principle, by running time backwards and using a shooting argument in a small neighbourhood of M_{1}, any heteroclinic orbit can be computed as accurate as the numerical time integrator allows. However we can be more precise and prepare the data even better by noticing that for p\equiv 0 the equation for q in (S) decouples completely from the rest and can be solved explicitly. Further, using this analytic value of q an exact solution can be also derived for r in the case when \displaystyle n=\frac{1}{k}, k\geq 1,\ k\in\mathbb{Z} :
\displaystyle\alpha=0,\quad p\equiv 0,\quad q(\eta)=\frac{1}{1+e^{\eta}},  (63)  
\displaystyle r(\eta)=\frac{r_{0}\left(1+e^{\eta}\right)^{k}}{\displaystyle% \sum_{j=0}^{k}\frac{kW_{0}}{kW_{0}j}\begin{pmatrix}k\\ j\end{pmatrix}e^{j\eta}},\quad\text{where}\ W_{0}=\frac{(m+n)r_{0}}{\lambda}. 
Step 2. (Projection step, \alpha=0,\ p\equiv 0) At this step we integrate numerically the equation for s using the exact values of q(\eta),r(\eta) found in the previous step. The integration can be performed by either AUTO or any other numerical integrator. The integration timespan is chosen to be \eta\in[\eta_{max},\eta_{max}] so that the starting point (p,q,r,s)_{\eta=\eta_{max}} and the ending point (p,q,r,s)_{\eta=\eta_{max}} both fall in small neighbourhoods of M_{0} and M_{1} respectively. In particular we choose the starting point so that (p,q,r,s)_{\eta=\eta_{max}}=M_{0}+\epsilon_{0}\nu_{0}, \nu_{0}=X_{02} and \epsilon_{0} as small parameter. Using this an initial value we integrate numerically the following nonautonomous equation for s
\dot{s}=s\Big{(}\frac{mn}{\lambda}(r(\eta)a)+q(\eta)\frac{1}{\lambda}r(% \eta)\big{(}s(1+m+n)\big{)}\frac{n}{\lambda}\Big{)}. 
At the end of the calculation we project the vector (p,q,r,s)_{\eta=\eta_{max}}M_{1} to the stable subspace of M_{1}. Indeed, we can find explicitly \epsilon_{1}\ll 1 and \nu_{1}\in\underset{}{\textrm{Span}}\{X_{11},X_{12},X_{14}\}, \nu_{1}=1 such that \pi\big{(}(p,q,r,s)_{\eta=\eta_{max}}M_{1}\big{)}=\epsilon_{1}\nu_{1}, where \pi denotes the projection. At the completion of Step 2 we have the solution (p,q,r,s)(\eta),\eta\in[\eta_{max},\eta_{max}] at discrete levels \eta_{i},i=0,\dots N for \alpha=0 and lying in the plane p\equiv 0.
Step 3. (Continuation with \alpha\neq 0,\ p=0) The goal in this step is to create a set of orbits in the plane p\equiv 0 but with \alpha not any more trivial. To that effect, we use the well prepared data obtained in Step 2 and we run AUTO with \nu_{0} fixed but allowing \alpha, m, \lambda, \epsilon_{0}, \epsilon_{1}, \nu_{1} be continued. The continuation process performed by AUTO creates a family of orbits with the following characteristics : a) emanate from a small neighbourhood of size \epsilon_{0} of M_{0} in the direction of \nu_{0}, b) terminate in a small neighbourhood of size \epsilon_{1} of M_{1}, c) lie in the plane p\equiv 0 but with \alpha\neq 0.
Step 4. (Continuation with \alpha\neq 0,\ p\neq 0) In this step we capture the desired heteroclinic orbit connecting M_{0} and M_{1}. From the family of orbits obtained in Step 3 we select one according to the physical relevance of the parameters \alpha, m, n and \lambda. We run AUTO again allowing \nu_{0} in \underset{}{\textrm{span}}\{X_{01},X_{02},X_{03}\} to be continued, thus leaving the plane p\equiv 0. AUTO generates a family of orbits emanating from M_{1}, terminating in a neighbourhood of size \epsilon_{0} of M_{0} and is the 2surface of heteroclinic orbits of W^{u}(M_{0})\cap W^{s}(M_{1}). One of these orbits is the desired heteroclinic orbit with \nu_{0}=X_{01}.
Remark 9.1.
We note that the exact solution of r in (63) is valid only for values of n of the form n=\frac{1}{k},\ k\geq 1,k\in\mathbb{Z}. In the case that n is not of this form then one can rely on an numerical integrator for solving as accurately as possible the equation for r using the exact value of q.
9.2 Numerical Results
In this section we illustrate the computation of the heteroclinic orbit of (S), following the steps described in detail above. Further, using the heteroclinic orbit (p,q,r,s) of system (S) we compute the associated selfsimilar solution in terms of the original variables v(x,t),\ u(x,t),\ \theta(x,t),\ \sigma(x,t).
We begin by giving the explicit relation of the variables (p(\eta),q(\eta),r(\eta),s(\eta)) to the original variables v(x,t),\ u(x,t), \theta(x,t),\ \sigma(x,t). Indeed, collecting the transformations and change of variables described in (24), (37) and (39), we obtain
\displaystyle v(x,t)  \displaystyle=\frac{1}{b}t^{b}\ \xi^{b_{1}}\ p^{\frac{(\alpha+m+n)}{D}}s^{% \frac{\alpha}{D}}r^{\frac{n}{D}}\quad\qquad u(x,t)=t^{b+\lambda}\ \xi^{b_{1}% 1}\ p^{\frac{(1+\alpha)}{D}}s^{\frac{\alpha}{D}}r^{1+\frac{n}{D}}  (64)  
\displaystyle\theta(x,t)  \displaystyle=t^{c}\ \xi^{c_{1}}\ p^{\frac{(1+m+n)}{D}}s^{\frac{m+n1}{D}}r^{% \frac{2n}{D}},\qquad\sigma(x,t)=t^{d}\ \xi^{d_{1}}\ p^{\frac{(\alpha+m+n)}{D% }}s^{\frac{\alpha}{D}}r^{\frac{n}{D}} 
where \eta=\log\xi,\ \xi=t^{\lambda}x and a,b,c,d,D as in (21), (22). We present now the results of three numerical experiments. All the computations where performed with \eta_{max}=10 and \lambda=\frac{1}{2}\lambda_{max} where \lambda_{max} is the upper bound of \lambda in (43). The initial values of \alpha, m are taken as \alpha=0 and m=0.6,\ 0.5,\ 0.5, respectively. The corresponding values of n remained fixed throughout the process and were n=0.025,\ 0.0125,\ 0.01. Following the process described in Section 9.1, the software AUTO was able to perform the continuation process and capture the desired heteroclinic orbit. The resulting values for \alpha and m are shown in the figures, along with the value of the parameter L_{p}=\alpha+m+n. The change of sign of L_{p} from positive to negative signals the onset of localization.
Figures 8, 9 and 10 illustrate the emergence of localization by depicting the profiles of the original variables v,\ u,\ \theta,\ \sigma at a few time instances. The vertical axes, except for the velocity v, are in logarithmic scale, however the corresponding yrange of values for each variable, is the same in all figures. In part (a) of each figure the velocity profile is depicted, which eventually, attains the shape of a step function. Parts (b) and (c) of the figures present the localization in strain rate and temperature respectively. In both cases the initial profile is a small perturbation of a constant state which at later time localizes at the origin. On the other hand, part (d) of figures shows the collapse of the stress to zero. The rate of localization at the origin differs for each numerical experiment and depends on the values of the materials parameters \alpha,\ m,\ n and L_{p}. In particular, the onset of localization is characterized by the parameter L_{p} taking a negative value. Further, the rate of localization is determined by the magnitude of this negative value, with larger negative values indicating faster localization, as it is observed in Figures 810.
Appendix A The loss of hyperbolicity for n=0
Consider the system (1) when n=0, that is the viscoplastic effects are neglected. Then (2) reads
\sigma=\tau(\theta,\gamma)=\theta^{\alpha}\gamma^{m}  (65) 
and (1) is written as a first order system
\begin{pmatrix}v_{t}\\ \theta_{t}\\ \gamma_{t}\end{pmatrix}=\underbrace{\begin{pmatrix}0&\tau_{\theta}(\theta,% \gamma)&\tau_{\gamma}(\theta,\gamma)\\ \tau(\theta,\gamma)&0&0\\ 1&0&0\\ \end{pmatrix}}_{\text{$\triangleq B(\theta,\gamma)$}}\begin{pmatrix}v_{x}\\ \theta_{x}\\ \gamma_{x}\end{pmatrix}.  (66) 
We check hyperbolicity for (66). The characteristic speeds are the roots of
\displaystyle\det\big{(}B\lambda\textrm{I}\big{)}  \displaystyle=\lambda\big{(}\lambda^{2}(\tau_{\theta}\tau+\tau_{\gamma})\big% {)}=0 
The system is thus hyperbolic when \tau_{\theta}\tau+\tau_{\gamma}>0 and elliptic in the tdirection when \tau_{\theta}\tau+\tau_{\gamma}<0. Observe that along the evolution of (66) and under the conditions for loading of interest in our problem, we have that \gamma is increasing; the equation
\theta_{t}=\tau(\theta,\gamma)\gamma_{t} 
implies that
\tau_{\theta}\tau+\tau_{\gamma}=\frac{d}{d\gamma}\tau(\theta,\gamma)\,. 
We conclude that the system is hyperbolic before the maximum of the stressstrain curve, and elliptic beyond the maximum point. For the constitutive law (65) a computation shows
\displaystyle\tau_{\theta}\tau+\tau_{\gamma}  \displaystyle=\theta^{\alpha}\gamma^{m1}\big{(}\alpha\frac{\gamma^{m+1}}{% \theta^{1+\alpha}}+m\big{)}  
\displaystyle=\frac{\alpha+m}{1+\alpha}+\frac{\alpha(1+m)}{\theta^{1+\alpha}}% \big{(}\frac{\theta_{0}(x)^{1+\alpha}}{1+\alpha}\frac{\gamma_{0}(x)^{1+m}}{1+% m}\big{)} 
In the region \alpha>m the stressstrain curve may be initially increasing (depending on the data) but eventually decreases.
The system (66) admits the class of uniform shearing solutions
v_{s}(x)=x\,,\quad\gamma_{s}(t)=t+\gamma_{0}\,,\quad\theta_{s}(t)\;\mbox{is % determined by solving}\;\begin{cases}\frac{d\theta_{s}}{dt}=\tau(\theta_{s},% \gamma_{s})&\\ \theta_{s}(0)=\theta_{0}&\\ \end{cases}  (67) 
where \gamma_{0}, \theta_{0}> are the initial strain and temperature, respectively. We linearize around the uniform shearing solution by setting
v=x+\hat{V}(t,x)\,,\quad\theta=\theta_{s}(t)+\hat{\Theta}(t,x)\,,\quad\gamma=% \gamma_{s}(t)+\hat{\Gamma}(t,x) 
and obtain the linearized system satisfied by the perturbation (\hat{V},\hat{\Theta},\hat{\Gamma}),
\begin{pmatrix}\hat{V}_{t}\\ \hat{\Theta}_{t}\\ \hat{\Gamma}_{t}\end{pmatrix}=B(\theta_{s}(t),\gamma_{s}(t))\begin{pmatrix}% \hat{V}_{x}\\ \hat{\Theta}_{x}\\ \hat{\Gamma}_{x}\end{pmatrix}+\begin{pmatrix}0&0&0\\ 0&\tau_{\theta}(\theta_{s},\gamma_{s})&\tau_{\gamma}(\theta_{s},\gamma_{s})\\ 0&0&0\\ \end{pmatrix}\;\begin{pmatrix}\hat{V}\\ \hat{\Theta}\\ \hat{\Gamma}\end{pmatrix}.  (68) 
The above calculation shows that, when \alpha>m, the linearized system loses hyperbolicity in finite time, past the maximum of the curve \sigma_{s}(t)t.
Appendix B The equilibria of the system (S)
We discussed in section 5 the equilibria M_{0} and M_{1} of (S). The remaining equilibria of (S) are listed below, and they are all functions of (\alpha,m,n,\lambda) that lie outside the the sector
\mathrsfs{P}=\{(p,q,r,s)\;\;p\geq 0,\,q\geq 0,\,r>0,\,s>0\} 
in the parameter range (23). The reader will find underlined the components indicating that the equilibrium lies outside the sector of interest: We recall the notations
a=\frac{2+2\alphan}{D}+\frac{2(1+\alpha)}{D}\lambda\,,\quad b=\frac{1+m}{D}+% \frac{1+m+n}{D}\lambda\,,\quad D=1+2\alphamn 
while t, t_{1} and t_{2} denote arbitrary real numbers.
\displaystyle\begin{array}[]{lllllll}(1)&\Big{(}0,&0,&\underline{0},&% \underline{0}&\Big{)},\\ (2)&\Big{(}0,&0,&\underline{0},&t&\Big{)}&\text{provided }\lambda=\frac{2% \alpha+2m+n}{2(\alphamn)},\\ (3)&\Big{(}0,&0,&\frac{n\alphaa(\alphamn)}{(1+\alpha)(m+n)},&\underline{0}&% \Big{)},\\ (4)&\Big{(}0,&1,&\underline{0},&\underline{0}&\Big{)},\\ (5)&\Big{(}0,&1,&\underline{0},&t&\Big{)}&\text{provided }\lambda=\frac{2% \alpha2mn}{1+m+n},\\ (6)&\Big{(}0,&1,&\frac{n\alphaa(\alphamn)}{(1+\alpha)(m+n)}+\frac{\lambda}{% m+n},&\underline{0}&\Big{)},\\ (7)&\Big{(}t,&0,&\underline{0},&\underline{0}&\Big{)}&\text{provided }\lambda=% \frac{2+2\alphan}{2(\alphamn)},\\ (8)&\Big{(}t,&1,&\underline{0},&\underline{0}&\Big{)}&\text{provided }\lambda=% \frac{22\alpha+n}{1+m+n},\\ (9)&\Big{(}t_{1},&0,&\underline{0},&t_{2}&\Big{)}&\text{provided $1+2\alpham% n=0$ and $\lambda=\frac{1m}{1+m+n}$},\\ (10)&\Big{(}t_{1},&1,&\underline{0},&t_{2}&\Big{)}&\text{provided $1+2\alpham% n=0$ and $\lambda=\frac{1m}{1+m+n}$},\\ \end{array} 
\displaystyle(11)\qquad\Big{(}  \displaystyle\underline{\tfrac{(\alphamn)(1+m+n)}{(1+\alpha)(1+m)}},\>% \tfrac{2(\alphamn)}{1+m}b,\>\tfrac{2(1+\alpha)}{1+2\alphamn},\>\tfrac{1+m+% n}{1+\alpha}\tfrac{n(1+2\alphamn)}{2(1+\alpha)^{2}}\Big{)},  
\displaystyle(12)\qquad\bigg{(}  \displaystyle\Big{(}\tfrac{2\alpha(1+m)}{D(1mn)}+\tfrac{2(\alphamn)}{D}% \lambda\Big{)}\Big{(}\tfrac{2\alpha(1+m)}{D(1mn)}\tfrac{1+m+n}{D}\lambda% \Big{)}\tfrac{1mn}{\lambda(2n)}\tfrac{1mn}{\lambda(1+m)},  
\displaystyle\qquad\left(\tfrac{2\alpha(1+m)}{D(1mn)}+\tfrac{2(\alphamn)}{% D}\lambda\right)\left(\tfrac{1+m}{D}+\tfrac{1+m+n}{D}\lambda\right)\tfrac{1m% n}{\lambda(1+m)},\ \tfrac{2n}{1mn},\ \underline{0}\bigg{)}. 
The generic equilibria in \mathrsfs{P} are M_{0}, M_{1}, (1), (34), (6), and (11 12); the rest are valid for specific parameter values.
Appendix C The linearized problems around M_{0} and M_{1}
The coefficient matrix for the linearized system (S) around the equlibrium M_{0} is
\displaystyle\begin{pmatrix}2&0&0&0\\ br_{0}&1&0&0\\ \frac{r_{0}}{n}(\lambda r_{0})&\frac{r_{0}}{n}&\frac{r_{0}}{n}\Big{(}\frac{% \alphamn}{\lambda(1+\alpha)}\frac{n\alpha}{\lambda(1+\alpha)r_{0}}\Big{)}&% \frac{r_{0}}{n}(\frac{\alpha r_{0}}{\lambda})\\ s_{0}(\lambda r_{0})&s_{0}&s_{0}\Big{(}\frac{\alphamn}{\lambda(1+\alpha)}+% \frac{n}{\lambda(1+\alpha)r_{0}}\Big{)}&s_{0}(\frac{r_{0}}{\lambda})\end{% pmatrix}=\begin{pmatrix}2&0&0&0\\ br_{0}&1&0&0\\ \frac{r_{0}}{n}(\lambda r_{0})&\frac{r_{0}}{n}&\frac{r_{0}}{n}\frac{1}{\lambda% }\Big{(}1s_{0}\frac{n}{r_{0}}\Big{)}&\frac{r_{0}}{n}(\frac{\alpha r_{0}}{% \lambda})\\ s_{0}(\lambda r_{0})&s_{0}&s_{0}\frac{1}{\lambda}(1s_{0})&s_{0}(\frac{r_{0}}% {\lambda})\end{pmatrix} 
The corresponding eigenvectors X_{0j} are collected in the matrix S_{0} as jth column vector, j=1,2,3,4.
\displaystyle S_{0}  \displaystyle=\begin{pmatrix}1&0&0&0\\ br_{0}&1&0&0\\ y_{1}&y_{2}&1&y_{4}\\ z_{1}&z_{2}&z_{3}&1\end{pmatrix},\quad\quad\begin{array}[]{l}\begin{pmatrix}y_% {1}\\ z_{1}\end{pmatrix}=(\lambda+b)r_{0}\begin{pmatrix}\frac{\frac{1+\alpha}{% \lambda}r_{0}+\frac{2}{s_{0}}}{\Delta_{1}}\\ \frac{\frac{n}{r_{0}}\big{(}\frac{1}{\lambda}+2\big{)}}{\Delta_{1}}\end{% pmatrix},\quad\begin{pmatrix}y_{2}\\ z_{2}\end{pmatrix}=\begin{pmatrix}\frac{\frac{1+\alpha}{\lambda}r_{0}+\frac{% \mu_{02}}{s_{0}}}{\Delta_{2}}\\ \frac{\frac{n}{r_{0}}\big{(}\frac{1}{\lambda}+\mu_{02}\big{)}}{\Delta_{2}}\end% {pmatrix}\\ z_{3}=n\bigg{(}\frac{\frac{1s_{0}}{\lambda}}{\frac{nr_{0}}{\lambda}+\frac{n% \mu_{0}^{+}}{s_{0}}}\bigg{)},\quad y_{4}=\frac{\frac{r_{0}}{\lambda}+\frac{\mu% _{0}^{}}{s_{0}}}{\frac{1s_{0}}{\lambda}},\end{array}  (69) 
where \Delta_{1}=\frac{1s_{0}}{\lambda}\big{(}\frac{1+\alpha}{\lambda}r_{0}+\frac{2% }{s_{0}}\big{)}\frac{n}{r_{0}}\big{(}\frac{1}{\lambda}+2\big{)}\big{(}\frac{r% _{0}}{\lambda}+\frac{2}{s_{0}}\big{)} and \Delta_{2}=\frac{1s_{0}}{\lambda}\big{(}\frac{1+\alpha}{\lambda}r_{0}+\frac{1% }{s_{0}}\big{)}\frac{n}{r_{0}}\big{(}\frac{1}{\lambda}+1\big{)}\big{(}\frac{r% _{0}}{\lambda}+\frac{1}{s_{0}}\big{)}. We find that y_{1},y_{2},y_{4}<0; z_{1},z_{2},z_{3}\sim{{O}}(n), provided n is sufficiently small.
Next, the coefficient matrix for the linearized system around M_{1} is
\displaystyle\begin{pmatrix}\frac{1+m+n}{\alphamn}&0&0&0\\ (b\lambda)r_{1}&1&0&0\\ \frac{r_{1}}{n}(\lambda r_{1})&\frac{r_{1}}{n}&\frac{r_{1}}{n}\Big{(}\frac{% \alphamn}{\lambda(1+\alpha)}\frac{n\alpha}{\lambda(1+\alpha)r_{1}}\Big{)}&% \frac{r_{1}}{n}(\frac{\alpha r_{1}}{\lambda})\\ s_{1}(\lambda r_{1})&s_{1}&s_{1}\Big{(}\frac{\alphamn}{\lambda(1+\alpha)}+% \frac{n}{\lambda(1+\alpha)r_{1}}\Big{)}&s_{1}(\frac{r_{1}}{\lambda})\end{% pmatrix}=\begin{pmatrix}\frac{1+m+n}{\alphamn}&0&0&0\\ (b\lambda)r_{1}&1&0&0\\ \frac{r_{1}}{n}(\lambda r_{1})&\frac{r_{1}}{n}&\frac{r_{1}}{n}\frac{1}{\lambda% }\Big{(}1s_{1}\frac{n}{r_{1}}\Big{)}&\frac{r_{1}}{n}(\frac{\alpha r_{1}}{% \lambda})\\ s_{1}(\lambda r_{1})&s_{1}&s_{1}\frac{1}{\lambda}(1s_{1})&s_{1}(\frac{r_{1}}% {\lambda})\end{pmatrix} 
In what follows we examine all possible cases: Except for the case \mu_{11}=\mu_{12}=1, four linearly independent eigenvectors are attained. In the exceptional case \mu_{11}=\mu_{12}=1 the repeated eigenvalue 1 has geometric multiplicity which is strictly less that its algebraic multiplicity.
As to the eigenvectors, notice that the eigenvalues for M_{1} (differently from those for M_{0}) have the chance to be repeated. The analysis below shows that, unless \mu_{11}=\mu_{12}=1, four linearly independent eigenvectors are attained. If the exceptional case takes place then we will supplement precisely one generalized eigenvector for the repeated eigenvalue 1.
Case 1. \frac{1+m+n}{\alphamn}\neq1; or \frac{1+m+n}{\alphamn}=1 but b=\lambda. This case yields four linearly independent eigenvectors. The eigenvectors X_{1j} are collected in the matrix S_{1} as jth column vector, j=1,2,3,4, and in the case of repeated eigenvalues the corresponding eigenvectors are understood as a basis for the associated subspace:
\displaystyle S_{1}  \displaystyle=\begin{pmatrix}1&0&0&0\\ x_{1}&1&0&0\\ y_{1}&y_{2}&1&y_{4}\\ z_{1}&z_{2}&z_{3}&1\end{pmatrix},\quad\quad\begin{array}[]{l}x_{1}=\begin{% cases}\frac{(b\lambda)r_{1}}{1+\mu_{11}}&\text{if $\mu_{11}\neq1$,}\\ 0&\text{otherwise,}\end{cases}\\ z_{3}=n\bigg{(}\frac{\frac{1s_{1}}{\lambda}}{\frac{nr_{1}}{\lambda}+\frac{n% \mu_{1}^{+}}{s_{1}}}\bigg{)},\quad y_{4}=\frac{\frac{r_{1}}{\lambda}+\frac{\mu% _{1}^{}}{s_{1}}}{\frac{1s_{1}}{\lambda}},\\ \end{array} 
\displaystyle\begin{pmatrix}y_{1}\\ z_{1}\end{pmatrix}=\begin{cases}(\lambda r_{1}+x_{1})\begin{pmatrix}\frac{% \lambda}{1s_{1}}\\ 0\end{pmatrix}&\text{if $\mu_{14}=\mu_{11}$,}\\ (\lambda r_{1}+x_{1})\begin{pmatrix}\frac{\frac{1+\alpha}{\lambda}r_{1}+\frac% {\mu_{11}}{s_{1}}}{\Delta_{3}}\\ \frac{\frac{n}{r_{1}}\big{(}\frac{1}{\lambda}+\mu_{11}\big{)}}{\Delta_{3}}\end% {pmatrix}&\text{otherwise,}\end{cases}\quad\begin{pmatrix}y_{2}\\ z_{2}\end{pmatrix}=\begin{cases}\begin{pmatrix}\frac{\lambda}{1s_{1}}\\ 0\end{pmatrix}&\text{if $\mu_{14}=\mu_{12}$,}\\ \begin{pmatrix}\frac{\frac{1+\alpha}{\lambda}r_{1}+\frac{\mu_{12}}{s_{1}}}{% \Delta_{4}}\\ \frac{\frac{n}{r_{1}}\big{(}\frac{1}{\lambda}+\mu_{12}\big{)}}{\Delta_{4}}\end% {pmatrix}&\text{otherwise,}\end{cases}  (70) 
where
\displaystyle\Delta_{3}  \displaystyle=\frac{1s_{1}}{\lambda}\big{(}\frac{1+\alpha}{\lambda}r_{1}+% \frac{\mu_{11}}{s_{1}}\big{)}\frac{n}{r_{1}}\big{(}\frac{1}{\lambda}+\mu_{11}% \big{)}\big{(}\frac{r_{1}}{\lambda}+\frac{\mu_{11}}{s_{1}}\big{)}  
\displaystyle=\frac{n}{r_{1}s_{1}}\det\left[\begin{pmatrix}\frac{r_{1}}{n}% \big{(}\frac{1s_{1}}{\lambda}\frac{n}{\lambda r_{1}}\big{)}&\frac{r_{1}}{n}% \frac{\alpha r_{1}}{\lambda}\\ s_{1}\frac{1s_{1}}{\lambda}&s_{1}\frac{r_{1}}{\lambda}\end{pmatrix}\mu_{11}% \textrm{I}\right]\neq 0,  
\displaystyle\Delta_{4}  \displaystyle=\frac{1s_{1}}{\lambda}\big{(}\frac{1+\alpha}{\lambda}r_{1}+% \frac{\mu_{12}}{s_{1}}\big{)}\frac{n}{r_{1}}\big{(}\frac{1}{\lambda}+\mu_{12}% \big{)}\big{(}\frac{r_{1}}{\lambda}+\frac{\mu_{12}}{s_{1}}\big{)}  
\displaystyle=\frac{n}{r_{1}s_{1}}\det\left[\begin{pmatrix}\frac{r_{1}}{n}% \big{(}\frac{1s_{1}}{\lambda}\frac{n}{\lambda r_{1}}\big{)}&\frac{r_{1}}{n}% \frac{\alpha r_{1}}{\lambda}\\ s_{1}\frac{1s_{1}}{\lambda}&s_{1}\frac{r_{1}}{\lambda}\end{pmatrix}\mu_{12}% \textrm{I}\right]\neq 0 
respectively for the corresponding cases.
Case 2. \frac{1+m+n}{\alphamn}=1 and b\neq\lambda: For this case \mu_{11}=\mu_{12}=1 has algebraic multiplicity two but its geometric multiplicity is one, so we replace the first column of S_{1} by the generalized eigenvector \big{(}\frac{1}{(b\lambda)r_{1}},0,y_{1}^{\prime},z_{1}^{\prime}\big{)}^{T}, where
\displaystyle\begin{pmatrix}y_{1}^{\prime}\\ z_{1}^{\prime}\end{pmatrix}=\begin{cases}\begin{pmatrix}\frac{\lambda}{1s_{1% }}\big{(}\frac{\lambda}{b\lambda}\frac{n}{r_{1}}z_{2}\big{)}\\ 0\end{pmatrix}&\text{if $\mu_{14}=1$,}\\ \frac{\lambda}{b\lambda}\begin{pmatrix}\frac{\frac{1+\alpha}{\lambda}r_{1}+% \frac{\mu_{11}}{s_{1}}}{\Delta_{3}}\\ \frac{\frac{n}{r_{1}}\big{(}\frac{1}{\lambda}+\mu_{11}\big{)}}{\Delta_{3}}\end% {pmatrix}+\frac{n}{r_{1}}\begin{pmatrix}\frac{y_{2}\big{(}\frac{r_{1}}{\lambda% }+\frac{\mu_{11}}{s_{1}}\big{)}+z_{2}\frac{\alpha r_{1}}{\lambda}}{\Delta_{3}}% \\ \frac{y_{2}\big{(}\frac{1s_{1}}{\lambda}\big{)}+z_{2}\big{(}\frac{1s_{1}}{% \lambda}+\frac{n}{r_{1}}\big{(}\frac{1}{\lambda}+\mu_{11}\big{)}\big{)}}{% \Delta_{3}}\end{pmatrix}&\text{otherwise.}\end{cases}  (71) 
\displaystyle 0  \displaystyle=Mat_{1}\begin{pmatrix}w\\ x\\ y\\ z\end{pmatrix}\mu\begin{pmatrix}w\\ x\\ y\\ z\end{pmatrix}=\begin{pmatrix}(\mu_{11}\mu)w\\ (b\lambda)r_{1}w+(\mu_{12}\mu)x\\ \frac{r_{1}}{n}\left[\lambda r_{1}w+x+\big{(}\frac{1s_{1}}{\lambda}\frac{n}{% r_{1}}\big{(}\frac{1}{\lambda}+\mu\big{)}\big{)}y+\frac{\alpha r_{1}}{\lambda}% z\right]\\ s_{1}\left[\lambda r_{1}w+x+\big{(}\frac{1s_{1}}{\lambda}\big{)}y\big{(}% \frac{r_{1}}{\lambda}+\frac{\mu}{s_{1}}\big{)}z\right]\end{pmatrix},  
\displaystyle A  \displaystyle\triangleq\begin{pmatrix}\frac{1s_{1}}{\lambda}\frac{n}{r_{1}}% \big{(}\frac{1}{\lambda}+\mu\big{)}&\frac{\alpha r_{1}}{\lambda}\\ \frac{1s_{1}}{\lambda}&\big{(}\frac{r_{1}}{\lambda}+\frac{\mu}{s_{1}}\big{)}% \end{pmatrix}\begin{pmatrix}y\\ z\end{pmatrix}=(\lambda r_{1}w+x)\begin{pmatrix}1\\ 1\end{pmatrix}  
\displaystyle A^{1}  \displaystyle=\frac{1}{\Delta}\begin{pmatrix}\big{(}\frac{r_{1}}{\lambda}+% \frac{\mu}{s_{1}}\big{)}&\frac{\alpha r_{1}}{\lambda}\\ \frac{1s_{1}}{\lambda}&\frac{1s_{1}}{\lambda}+\frac{n}{r_{1}}\big{(}\frac{1% }{\lambda}+\mu\big{)}\end{pmatrix},\quad\Delta=\frac{1s_{1}}{\lambda}\big{(}% \frac{1+\alpha}{\lambda}r_{1}+\frac{\mu}{s_{1}}\big{)}\frac{n}{r_{1}}\big{(}% \frac{1}{\lambda}+\mu\big{)}\big{(}\frac{r_{1}}{\lambda}+\frac{\mu}{s_{1}}\big% {)} 
Acknowledgement. The authors thank Prof. Peter Szmolyan for valuable discussions on the use of geometric singular perturbation theory.
References
 [1] M. Bertsch, L. Peletier, and S. Verduyn Lunel, The effect of temperature dependent viscosity on shear flow of incompressible fluids, SIAM J. Math. Anal. 22 (1991), 328–343.
 [2] R.J. Clifton, High strain rate behaviour of metals, Appl. Mech. Rev. 43 (1990), S9S22.
 [3] R.J.Clifton, J.Duffy, K.A.Hartley and T.G.Shawki, On critical conditions for shear band formation at high strain rates, Scripta Met., 18 (1984), pp. 443448.
 [4] L.S.Costin, E.E.Crisman, R.H.Hawley, and J.Duffy, On the localization of plastic flow in mild steel tubes under dynamic torsional loading In ”Proc. 2nd Conf. on the Mechanical Properties of Materials at high rates of strain”, Inst. Phys. Conf. Ser. no 47, Oxford, 90, 1979.
 [5] C.M. Dafermos and L. Hsiao, Adiabatic shearing of incompressible fluids with temperaturedependent viscosity. Quart. Applied Math. 41 (1983), 45–58.
 [6] E.J. Doedel, AUTO: A program for the automatic bifurcation analysis of autonomous systems, Cong. Numer. 30 (1981), 265–284.
 [7] E.J. Doedel and J.P. Kernevez, AUTO: Software for continuation and bifurcation problems in ordinary differential equations, Applied Mathematics Report, California Institute of Technology, (1986), 226 pages.
 [8] E.J. Doedel, A.R. Champneys, T.F. Fairgrieve, Y.A. Kuznetsov, B. Sandstede, and X. Wang, AUTO 97: Continuation And Bifurcation Software For Ordinary Differential Equations (with HomCont), (1999), (http://indy.cs.concordia.ca/auto/)
 [9] N. Fenichel, Persistence and smoothness of invariant manifolds for flows, Indiana Univ. Math. J. 21 (1972) 193–226.
 [10] N. Fenichel, Asymptotic stability with rate conditions, Indiana Univ. Math. J. 23 (1974) 1109–1137.
 [11] N. Fenichel, Asymptotic stability with rate conditions II, Indiana Univ. Math. J. 26 (1977) 81–93.
 [12] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, J. Differ. Equations 31 (1979), 53–98.
 [13] C. Fressengeas and A. Molinari Instability and localization of plastic flow in shear at high strain rates, J Mech Phys Solids 35 (1987), pp. 185211.
 [14] K.A.Hartley, J.Duffy, and R.J.Hawley, Measurement of the temperature profile during shear band formation in steels deforming at highstrain rates J. Mech. Physics Solids, 35 (1987), 283301.
 [15] J.W. Hutchinson and K.W. Neale, Influence of strainrate sensitivity on necking under uniaxial tension, Acta Metallurgica 25 (1977), 839846.
 [16] Th. Katsaounis, J. Olivier, and A.E. Tzavaras, Emergence of coherent localized structures in shear deformations of temperature dependent fluids, Archive for Rational Mechanics and Analysis 224 (2017), 173–208.
 [17] Th. Katsaounis and A.E. Tzavaras, Effective equations for localization and shear band formation, SIAM J. Appl. Math. 69 (2009), 1618–1643.
 [18] Th. Katsaounis, M.G. Lee, and A.E. Tzavaras, Localization in inelastic rate dependent shearing deformations, J. Mech. Phys. of Solids 98 (2017), 106–125.
 [19] M.G. Lee, Th. Katsaounis, and A.E. Tzavaras, Localization of Adiabatic Deformations in Thermoviscoplastic Materials, In Proceedings of the 16th International Conference on Hyperbolic Problems: Theory, Numerics, Applications (HYP2016), to appear.
 [20] C. K. R. T. Jones, Geometric singular perturbation theory. Dynamical systems (Montecatini Terme, 1994), pp 44Ã118, Lecture Notes in Math., 1609, Springer, Berlin, 1995.
 [21] C. Kuehn, Multiple time scale dynamics, Applied Mathematical Sciences, Vol. 191 (Springer Basel 2015).
 [22] M.G. Lee and A.E. Tzavaras, Existence of localizing solutions in plasticity via the geometric singular perturbation theory, Siam J. Appl. Dyn. Systems 16 (2017), 337–360.
 [23] A.Molinari and R.J.Clifton, Analytical characterization of shear localization in thermoviscoplastic materials J. Appl. Mech., 54 (1987), 806812.
 [24] T.G. Shawki and R.J. Clifton, Shear band formation in thermal viscoplastic materials, Mech. Mater. 8 (1989), 13–43.
 [25] P. Szmolyan, Transversal heteroclinic and homoclinic orbits in singular perturbation problems, J. Differ. Equations 92 (1991), 252–281.
 [26] A.E. Tzavaras, Shearing of materials exhibiting thermal softening or temperature dependent viscosity, Quart. Applied Math. 44 (1986), 1–12.
 [27] A.E. Tzavaras, Effect of thermal softening in shearing of strainrate dependent materials. Archive for Rational Mechanics and Analysis, 99 (1987), 349–374.
 [28] A.E. Tzavaras, Plastic shearing of materials exhibiting strain hardening or strain softening, Arch. Ration. Mech. Anal. 94 (1986), 39–58.
 [29] A.E. Tzavaras, Nonlinear analysis techniques for shear band formation at high strainrates, Appl. Mech. Rev. 45 (1992), S82–S94.
 [30] T.W. Wright, The Physics and Mathematics of Shear Bands. (Cambridge Univ. Press 2002).
 [31] C. Zener and J. H. Hollomon, Effect of strain rate upon plastic flow of steel, J. Appl. Phys. 15 (1944), 22–32.