Soliton solutions in geometrically nonlinear Cosserat micropolar elasticity with large deformations

We study the fully nonlinear dynamical Cosserat micropolar elasticity problem in space with three dimensionals with various energy functionals dependent on the microrotation $\overline{R}$ and the deformation gradient tensor $F$ . We derive a set of coupled nonlinear equations of motion from first principles by varying the complete energy functional. We obtain a double sine-Gordon equation and construct soliton solutions. We show how the solutions can determine the overall deformational behaviour and discuss the relations between wave numbers and wave velocities thereby identifying parameter values where the waves cannot propagate.


Introduction
Classical elasticity is based on considering materials whose idealised material points are structureless. Any possible internal properties are neglected in the classical theory. A microcontinuum, on the other hand, is a continuous collection of deformable materials points [1,2,3,4,5]. The characteristic aspect of the theory with microstructure is that we assume the microelements to exhibit an inner structure attached to so-called directors, which span an internal three-dimensional space. These can, for instance, rotate and deform. The most general case of this microcontinuum is the micromorphic continuum which has nine additional degrees of freedom when compared with classical elasticity theory. These additional degrees of freedom consist of 3 microrotations, 1 (micro) volume expansion and 5 (micro) shear deformations of the directors.
A special model arises when the extra degrees of freedom are reduced to rigid rotations. This means only 3 additional microrotations are considered, in addition to the classical translational deformation field. This theory is often referred to as Cosserat elasticity or micropolar elasticity and was originally proposed in full generality by the Cosserat brothers in 1909, see [6]. In some ways, their work was ahead of their time and was consequently largely forgotten for many decades. Starting from the 1950s interests in this theory increased and many advances were made since then [7,8,9,10,11,12,13,14,15,16].
We denote the microrotation by R = exp(X), where X is a one-parameter subgroup generated by the angle φ = φ(x, t) and the deformation gradient vector is related to the displacement vector u as F = ∇ϕ = 1+∇u. The deformation gradient F can be written using the polar decomposition F = RU and we can express u as a function of ψ = ψ(x, t), see Fig. 1 Once one can identify and collect the relevant energy functionals for the Cosserat micropolar elasticity, the equations of motion for the system can be found through the corresponding Euler-Lagrange equation. Due to the highly nonlinear nature of the system, various attempts were made to simplify the process under relatively weak restrictions and a simple ansatz. Spinor methods were used in [17] to simplify the Euler-Lagrange equation and subsequent works appeared in [18,19], with an intrinsically two-dimensional model studied in [20]. An investigation in optimisation of the Cosserat shear-stretch energy in searching for the optimal Cosserat rotation is made in [21,22]. In [23], the polarity of ferromagnets gave rise to the description of the defects in order parameters as the solitary waves under the external magnetic stimuli, followed by the study in the elastic crystals as a micropolar continuum in [24], again with the description of soliton solution for the topological defects. The reader may imagine an artificial discrete model where rotational discs connected by line and torsional springs are coupled in chiral fashion as a mechanical realisation of the model. Variants of the geometrically nonlinear Cosserat model are also used to describe lattice rotations in metal plasticity, see for instance [25].
In the recent paper [26], the dynamical Cosserat model was investigated by analysing the geometrically nonlinear and coupled nature of the system, in which the linearised energy functionals are used to simplify the problem significantly. It allowed the reduction of the coupled system of PDEs to a sine-Gordon equation, which in turn yielded a soliton-like solutions both in rotational and displacement deformations under the assumption that displacements are small while large and multiple rotations are allowed.
In this paper, we present the solutions of elastic and rotational propagation of deformations in the complete dynamical Cosserat problem. This involves the total energy functional given by (1.1) We will start with exactly the same ansatz used in [26] such that the displacement deformation wave is a plane wave in the form of ψ = g(z − vt) for some arbitrary function g with wave speed v. We expect to obtain a similar system of equations but with additional nonlinear terms. It turns out that these will yield a double sine-Gordon type equation.
The primary mathematical interests in finding the equations of motion using the variational calculus come from the fact that many terms in the energy functional contain quantities such as R T Curl R, R T polar(F ), or R T F throughout the calculations. Since in general the elements R ∈ SO(3) do not commute, their variations require a careful treatment in the calculations. The plan of the paper is the following. In Section 2, after stating each energy functional in terms of F and R, we vary the total energy functional including kinetic energy. We collect terms from the variational field expressions with respect to F and R in Section 3 to obtain the complete coupled system of equations of motion. It turns out that if we impose the previous restriction, i.e. small displacements, the newly generated nonlinear coupling terms in the complete description are indeed responsible for the contribution in the additional terms of the sine-Gordon type equation, as shown in Section 4. This observation reduces the problem to solving the so-called double sine-Gordon equation [27] of a single function of φ = φ(z, t). In the final Section we illustrate the effects of rotational and displacement propagations in the simple model of microcontinuum with additional features of kink-antikink form of solutions and the profiles of the wave number k and wave velocity v relations.
rotation matrix, microrotation X skew-symmetric matrix generating R

The complete dynamical Cosserat problem
We introduce each energy functional for the full treatment of the geometrically nonlinear Cosserat problem in three-dimensional space. We will subtract the kinetic energies from relevant energy functionals before deriving the equations of motion. First, the energy functional for elastic deformations is where λ and µ are the standard Lamé parameters. The microrotations are governed by the energy functional V curvature defined by where κ i are the elastic constants for the microrotations. An interaction between elastic displacements and microrotations is described by the irreducible parts of the elastic deformations and microrotations, such as R T F − 1 and R T Curl R respectively to form the energy functional V interaction (F, R) defined by where χ 1 and χ 3 are the coupling constants. Finally, we will consider the Cosserat coupling term which is given by where µ c is the Cosserat couple modulus. The variations of the complete energy functional are quite involved. All required results are stated explicitly in Appendix A. Gathering all the variational terms (A.5), (A.13), (A.15) and (A.19), we will obtain the complete variational functional of the theory for the dynamical case The next step will be collecting the various expressions with respect to F and R to construct the field equations.

Displacements and rotations in one axis
Let us assume that the points in our continuum can only experience rotations about one axis, say the z-axis, which means we can choose The variation of this is simply In principle, the rotational and elastic waves can be either longitudinal or transverse in each case, hence four different combinations are possible. Here, we consider solutions in which both waves are longitudinal about the same axis, the z-axis in this case, so that we can write ψ = ψ(t, z) and φ = φ(t, z).
Further, we collect the relevant terms with respect to F and R separately. Unlike the case of δR, in which the variational kinetic term is readily written with respect to R, the variational kinetic term from the interaction energy functional is written with respect to δu. But the variation with respect to F can be restated as the variation with respect to u, hence with respect to ψ as we will see shortly.
Collecting terms for δF from (2.6) gives Now, the terms which appear in the variation with respect to F can be transformed into the variation with respect to ∇u, for any matrix A, as follow.
up to a boundary term. In this case, the contribution comes only from A 33 and we obtain We now include the kinetic variational term ρü δu = ρ∂ tt ψ δψ to obtain the equation of motion In the same way, we collect terms for δR to obtain Applying B : δR gives which is (3.11) Therefore, from (3.8) and (3.11), we obtain two equations of motion by varying the total energy functional with respect to F and R, respectively, as follows These can be written in component form as From this, we can see immediately that we will recover the result obtained in [26] if we assume the linearised energy functionals which lead to the approximations such as λφ 1 and µφ 1, while the matrix elements M remain unchanged.
The revised results of [23] were stated in [1], in which case the longitudinal wave is expressed as U (x, t) along the x axis with the rotational deformation φ(x, t) about x axis. The equations of motion are described as a system of coupled expressions, (3.15) where α, λ , µ , κ are isotropic material moduli used in [1] and Since the matrix N is diagonal, we do not have second order coupling terms in the equations of motion. And under the small displacement limit, the system is readily solvable using the conventional method for the one-dimensional d'Alembert's solution subject to the appropriate boundary conditions.

Solution for the double sine-Gordon equation
We assume that the elastic and rotational waves propagate with the same wave speed v and ψ = g(z − vt), so that ψ satisfies ∂ tt ψ = v 2 ∂ zz ψ. Without this assumption we are not able to construct a solution. Now, we define v 2 rot = M 11 and v 2 elas = M 22 . Then (3.12b) becomes Integrating with respect to z once gives in which we set the constant of integration to zero by imposing the boundary condition ψ = ∂ z ψ = 0 as z → ±∞. Substituting (3.17) and (3.18) into the remaining equation of motion (3.12a) gives Moreover, if we rescale z as then (3.19) reduces to, the so-called double sine-Gordon equation where The apparent singularity in b as v 2 approaches v 2 elas can be removed if we make the further trans- We note that this transformation on v would not change our assumption on ψ along with the rescaling on z, since ∂ tt ψ = v 2 ∂ zz ψ implies ∂ tt ψ =v 2 ∂ẑẑψ. The general solution of (3.21) is given in [27] as in which u must satisfy two conditions The simplest solution is of the form with Now, we can write the solution φ using the identity arcsin(x) = 2 acrtan where φ 0 is the rotational propagation solution based on the linearised energy functionals with corresponding k 0 and m 0 given by Now, consider the quantity for δ = ln 1 2 . We would like to see if this agrees with the argument of the exponential in (3.29). This can be done if we apply the reverse rescaling (3.20) of z and inverse transformation (3.23) of v. After some calculations, we obtain For the current case, by following the same reasoning we find that the rescaled variables and original variables are interchangeable by the expression with δ = ln 1 2 . We must notice that the matrix M used in (3.35) and (3.30) is the same (3.14). The Lamé parameters λ and µ are brought into play in the fully nonlinear case through the quantity m 2 , while those parameters are missing in m 2 0 when considering the approximations λφ 1 and µφ 1. Consequently, we have to treat a more complicated form of k with an additional contribution from b. And it is clear that we can recover the solution (3.33) if we apply the restrictions λφ 1 and µφ 1, which will effectively lead to b = 0 and m → m 0 . For ψ, first we write X (hence u) in terms of z and v.
Plugging (3.36) into (3.17) gives, If we put s = z − vt, then this becomes a second-order ordinary differential equation for g(s).
We integrate twice with respect to s using the boundary conditions ψ (±∞, t) = ψ(±∞, t) = 0 to obtain (3.40) The constant C is Using the restriction λφ 1 and µφ 1, these solutions reduce to the one we obtained in [26] φ 0 = 4 arctan e ±k0(z−vt)±δ (3.42) In Fig. 3, the soliton solutions for φ(z, t) and ψ(z, t) are given at t = 0 with corresponding values of k. As the rotational wave φ(z, t) propagates with a speed v along the z-axis, the points of microcontinuum (displayed as pendulums along the z-axis) experience microrotational deformations perpendicular to the axis. In the same way the longitudinal solution ψ(z, t) gives rise to the compressional deformation wave propagating with the same speed v, on the points of macrocontinuum (shown as beads) along the axis. As we vary the values of k, the widths of the soliton solutions are changed and this affects the overall deformational behaviours both in rotation and displacement.

Properties of solutions
We notice that there might be possible singularity issues in the amplitude of ψ(z, t) in (3.39) as v 2 approaches v 2 elas . In order to resolve this problem, we would like to look closely at k as a function of v taking account of all nine parameters, {κ 1 , κ 3 , χ 1 , χ 3 , ρ, ρ rot , µ c , λ, µ}. We consider only the positive roots of k 2 to understand the possible range of k for a given v. After putting all relevant parameters in (3.35), we obtain Now, to determine whether k possesses any singularity, we compute the discriminant of the quartic of v in the denominator regarding it as a quadratic equation for v 2 . ( where we put v 2 χ ≡ M 12 . This is strictly non-negative, so that we can have four roots of v in the denominator of (4.1), which will cause the singularity of k. We denote the four distinct roots as v i , i = 1, 2, 3, 4 and assume that v 1 < v 2 < 0 < v 3 < v 4 . In particular, we write explicitly The Now, we plot the profiles of v as a function of k, this is given implicitly by (4.1), and we consider only the positive values of v for the simplicity. At this time, we only have two asymptotic lines of v 3 and v 4 (again we assume v 3 < v 4 ). And we assume that v elas > v rot .
Two characteristic types of parameter ranges for v with various values for a set of parameters with relevant asymptotic lines and the locations of v 0 , v elas and v rot are given in Fig. 4. The dominating set of parameters in determining the characteristics is the set of constants {λ, µ, µ c } of the energy functional V elastic . Notably, we observe that we only alter the value of the parameter µ c to obtain the type (b) solution from the type (a) solution while keeping all remaining parameters unchanged. The values of v elas and v rot are located inside (or on the boundary of) the shaded region surrounded by asymptotic lines, which can be shown directly from (4.3). The threshold in transition from the type (a) to (b) is evidently the relative positions between v 0 and v 4 . If v 0 > v 4 we will have the type (a) and if v 0 < v 4 then the type (b).
In both types (a) and (b) solutions, there exist regions (the shaded regions) in which v cannot be defined for a given k, solutions with such parameter choices do not exist. In case of type (a), the values of v are defined in v ∈ [0, v 3 ) and v ∈ (v 4 , v 0 ]. The upper limit of v is bounded by v 0 and we can see that v 0 → ∞ as µ c → 0 which is evident from (4.4), see Fig. 5.
On the other hand, for the type (b), the position of v 0 is v 3 < v elas < v 0 < v 4 . Now, the line of v 0 acts the role of the boundary line along with v 3 in (b). So v takes the values in the region v ∈ [0, v 3 ) and v ∈ [v 0 , v 4 ). We must notice that for type (b) solutions, the value of v 0 cannot be exactly v elas due to the restriction (4.4), as long as we have nonzero λ. We observe that v 0 approaches v elas as µ c → ∞, but the lower profile of v in (b) will be shifted to the right indefinitely, i.e. k → ∞, see Fig. 5. In the limit µ c → ∞, it is clear that we will have a profile of type (b). Also we can see from (3.35) that m 2 → m 2 0 , hence k 2 → k 2 0 . This suggests that b becomes negligible and we will be left with the soliton solution φ → φ 0 of the form (3.29).
Next, we consider the limit ρ rot ρ v 4 1.  Fig. 4, is pushed up to the infinity (left). In the limit µ c → ∞ the lower profile of type (b) will shift to infinity along the k axis (right).
In this limit, we can approximate the expressions of v 3 and v 4 given by (4.3) as follows Hence we can see that v elas approaches to v 4 and v rot approaches to v 3 for the type (a) parameter choice. In case of type (c) of Fig. 6, we set v χ = 0 (i.e., 3χ 1 − χ 3 = 0) to illustrate that v elas = v 4 and v rot = v 3 and that the lines of v elas and v rot play the role of asymptotic lines. In this case, the matrix M of (3.14) becomes diagonal and the system looks similar to (3.15). Of course, if we had assumed that v elas < v rot , then we would have v rot = v 4 and v elas = v 3 . We may obtain the similar observation in type (b) diagram by adjusting µ c , but v elas , v rot → v 3 . Furthermore, in the same limit of v χ = 0, if we set an additional condition that v elas = v rot , then we will have one asymptotic line v elas as shown in the diagram, type (d) of Fig. 6 and the matrix M will simply become the identity matrix (up to the rescaling). Now, the amplitude of ψ in (3.39) is determined by two coefficients (the matrix element M 21 can be written in terms of v 2 χ ≡ M 12 ), . The analytic investigation on the profiles of v as a function of k provides us the clue that the amplitude of ψ cannot be arbitrarily large. As k → ∞, we have v 2 → v 2 elas but the statement that the value of v 2 elas approaches v 2 4 is equivalent to say that v 2 χ → 0, as we can see directly from (4.7). Hence the first coefficient in (4.7) is assumed to remain finite in this limit. Similarly, the second coefficient cannot be arbitrarily large. For given k and (v 2 − v 2 elas ) will compensate each other as k → ∞. This is shown in the type (c), or more extreme case, the type (d) in Fig. 6.

Conclusion
We extended the previous study of the deformations considered in [26] to include the fully nonlinear model with arbitrarily large rotations and displacements. This discussion gave us further insights into the nature of the nonlinear geometry of Cosserat micropolar elasticity. The solution φ differs from φ 0 via the different form of k in (3.35). On the other hand, the displacements ψ and ψ 0 differ by additional nonlinear terms.
The soliton solutions for both rotations and displacements were obtained from the equations of motion and these allowed us to understand the geometric interpretation of the deformation waves. The physically dominant parameters of the complete model are the Lamé parameters {λ, µ} and the Cosserat couple modulus µ c . This becomes evident by looking at the k dependency, or equivalently m dependency, of the soliton solutions on these parameters.
The various values for k in the soliton solutions for φ and ψ give different overall behaviours while other values of parameters are fixed. Regarding the microrotations, the effect becomes apparent for large values of k, which induce high-frequency of localised energy distribution on the narrow width affected cross section both for the rotational and displacement deformations, whereas small values of k induce gradual and broad energy distribution for the deformations over the microcontinuum media. The role of k can be understood using a simple model of beads and pendulums as shown in Fig. 3.
A consideration for the deformation waves of higher dimensions would be a natural extension of the procedure. Some other candidates for further applications would include an investigation of domain walls in topological defects (e.g. ferromagnets) in connection with micropolar deformation. Moreover vortices as topological solitons with a notion of spontaneous symmetry breaking as a phase transition by Cosserat elasticity would be also be an interesting subject of study.

A Variations of energy functional
We would like to vary each energy functional using some of identities listed in Notation and Appendix. First, for V elastic , we can expand the expression using the definition of X 2 = X, X = tr(XX T ) and symM = 1/2(M + M T ) as Variation of this is If we want to study the dynamical problem, we must take the kinetic term into account in the elastic energy functional.
where ρ is the constant density and ϕ is the deformation vector. If we vary this term we will obtain δV elastic,kinetic = −ρφ δϕ.
But, since ∇ϕ = 1 + ∇u implies δϕ = δu andφ =ü, the variation of elastic kinetic term can be rewritten as δV elastic,kinetic = −ρü δu and the variation of dynamical expression for the elastic energy functional becomes Similarly, for the curvature functional, we can expand it as This is a functional dependent only on R, but the actual variation will involve rather complicated quantites such as δ Curl R multiplied by a tensor. To overcome this problem, we introduce the following identity. Let A(R) and B(R) be two matrix valued functions depending on the rotation R. Then, by direct calculation, one can show that an identity for any rank-two tensors A and B, In this way, we find the variation of curvature term δV curvature (R) = (κ 1 −κ 2 ) (Curl R)R T (Curl(R))+Curl R(Curl R) T R +(κ 1 +κ 2 ) Curl Curl R Again, for the dynamical case, we need to include the kinetic term defined as V curvature,kinetic = ρ rot Ṙ 2 = ρ rot tr(ṘṘ T ) (A. 12) with variational form given by δV curvature,kinetic = −2ρ rotR : δR. Therefore, the variation of dynamical expression for the curvature energy functional can be written as δV curvature (R) = (κ 1 −κ 2 ) (Curl R)R T (Curl(R))+Curl R(Curl R) T R +(κ 1 +κ 2 ) Curl Curl R For the interaction energy functional, we expand terms dev sym(R T Curl R) and dev sym(R T F − 1) to write 14) The variation of this involves the quantity δ Curl R as in the case of V curvature , so we use the identity (A.9) to obtain We note that this depends on R and R = polar(F ), hence depends on R and F . Therefore, the variation of coupling energy functional is of the form The term in the brackets in the second term can be written as In the first step, we used the chain rule and in the second and last steps we used the identities given in Appendix. Then the variation of coupling energy becomes